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I. 


INTRODUCTION 


Frequency  domain  beamforming  is  accomplished  by  applying 
appropriate  phase  shifts  at  the  sensor  outputs  of  an  array 
to  account  for  the  relative  propagation  delays  of  a  signal 
from  a  particular  direction.  The  phase-shifted  signals  from 
all  sensors  are  then  added  together  coherently  to  realize 
the  full  array  gain.  Discrete  Fourier  Transform  (DFT)  beam¬ 
forming  is  the  usual  method  of  determining  the  direction  of 
arrival  of  a  plane  wave  signal.  A  discrete  number  of 
spat-ial  frequency  bins  are  formed  and  each  bin  corresponds 
to  a  discrete  direction.  If  the  number  of  spatial  frequency 
bins  is  large,  very  fine  spatial  resolution  can  be  obtained. 

The  phase  shifts  needed  to  cancel  the  relative  propaga¬ 
tion  delays  can  be  determined  adaptively.  The  complex  LMS 
adaptive  algorithm  is  used  in  this  thesis.  The  LMS  adaptive 
filter  adjusts  its  adaptive  weights  recursively  to  minimize 
the  mean  square  difference  between  a  reference  signal  and 
its  estimate.  When  the  beam  is  steered  toward  a  signal 
propagating  in  a  particular  direction,  the  phase  of  the 
signals  at  all  sensors  must  be  the  same.  Therefore,  the 
signal  at  any  sensor  can  be  used  as  a  reference  which  the 
others  will  be  matched.  The  estimated  signal  is  obtained  by 
weighting  the  input  signal  by  the  current  adaptive  weights. 
Note  that  no  prior  knowledge  of  the  reference  signal  is 


required.  The  response  of  the  LMS  adaptive  filter  converges 
to  the  discrete  Wiener  filter  without  a  priori  knowledge  of 
the  input  [Refs.  1,2] .  [Ref.  1]  proposed  the  complex  LMS 
algorithm  to  deal  with  complex  inputs.  [Ref.  5]  addressed 
the  implementation  of  the  LMS  adaptive  filter  in  the  fre¬ 
quency  domain.  [Ref.  4]  used  the  LMS  adaptive  filter  in  the 
frequency  domain  to  estimate  the  bearing  of  a  plane  wave 
due  to  an  acoustic  source  radiating  a  sinusoidal  signal.  In 
this  application,  the  LMS  adaptive  filter  was  implemented  to 
estimate  the  phase  difference  between  two  sonar  arrays 
separated  by  a  distance  many  times  the  signal  wavelength. 

The  angle  of  arrival  of  a  plane  wave  can  be  estimated  if  the 
frequency  of  the  acoustic  signal  and  the  speed  of  wave  propa¬ 
gation  are  known  or  can  be  extracted  from  the  received  signal 
itself . 

The  objective  of  this  thesis  is  to  extend  the  results  in 
[Ref.  3]  and  [Ref.  4]  to  a  planar  array  of  M  *N  elements 
(hydrophones)  where  M  and  N  are  greater  than  two.  Such  an 
array  has  an  overall  size  many  times  the  wavelength  of  the 
received  signal.  The  inter-element  separation,  however,  is 
usually  maintained  at  a  distance  of  less  than  or  equal  to 
one-half  of  the  expected  minimum  wavelength.  This  requirement 
[Ref.  5]  prevents  the  occurrence  of  grating  lobes  in  the 
far-field  beam  pattern.  A  two-dimensional  array  allows 
spatial  resolution  in  both  azimuth  and  elevation.  Even 
though  the  detection  range  in  underwater  acoustics  is  large 


compared  to  the  ocean  depth,  the  effect  of  ray  bending  due 
to  the  inhomogeneous  ocean  medium  can  bend  the  incident 
acoustic  wave  such  that  it  can  appear  to  arrive  at  a  steeper 
or  shallower  angle  than  the  line  of  sight  angle  in  the  homo¬ 
geneous  medium  case.  The  elevation/depression  angle  is 
at  present  used  to  estimate  Convergence  Zone  (CZ)  and  Bottom 
Bounce  (BB)  ranges. 

In  this  thesis,  the  problem  of  a  plane  wave  incident  upon 
a  planar  array  of  M  * N  elements  is  studied.  The  acoustic 
wave  signal  at  each  of  the  elements  in  the  array  are  identical 
if  the  array  is  steered  in  the  direction  of  the  incident 
wavefront.  However,  if  the  main  lobe  of  the  array  is  not 
steered  properly,  the  plane  wave  signal  will  have  the  same 
spectral  content  at  each  element  but  modified  by  a  phase 
shift  proportional  to  the  location  of  the  element  with  respect 
to  some  reference  element.  These  undesirable  phase  shifts 
can  be  cancelled  by  applying  appropriate  phase  weights  at 
each  element  and  thereby  cophasing  the  total  array  output  to 
realize  its  array  gain.  [Ref.  4]  demonstrated  that  the  LMS 
adaptive  filter  can  achieve  phase  alignment  between  a  refer¬ 
ence  signal  and  input  signal  in  the  frequency  domain  by 
direct  application  of  the  complex  LMS  algorithm.  This  is 
equivalent  to  a  tapped  delay  line  structure  in  the  time  do¬ 
main.  However,  in  the  frequency  domain,  a  time  delay  t 
corresponds  to  multiplication  of  a  complex  number  that  is 
equal  to  e^wT  where  w  is  the  signal  frequency.  The 


implementation  of  the  LMS  adaptive  filter  in  the  frequency 
domain  requires  fewer  computations  per  iteration  than  in  the 
time  domain.  An  added  advantage  of  using  phase  weighting 
is  that  a  continuous  range  of  spatial  directions  can  be 
described,  whereas  in  a  tapped  delay  line  structure,  only 
finite  increments  of  delays  can  be  applied.  Figure  1  shows 
a  functional  diagram  of  an  N-element  adaptive  array  imple¬ 
mented  in  the  frequency  domain  [Ref.  6] . 

Chapter  II  of  this  thesis  describes  the  specific  struc¬ 
ture  of  the  adaptive  filter  and  the  equations  implemented  for 
simulations.  The  assumptions  made  in  the  model  are  discussed 
and  justified.  Several  modifications  to  the  complex  LMS 
adaptive  algorithms  are  made  to  increase  the  array's  spatial 
coverage  and  to  ensure  that  the  steady  state  phase  weights 
do  correspond  to  the  direction  of  the  incident  wave. 

In  Chapter  III,  a  passive  sonar  system  is  modeled  to 
test  the  ability  of  the  modified  complex  LMS  algorithm  to 
estimate  the  direction  of  a  source  in  the  presence  of  noise. 

The  simulation  program  is  implemented  in  VS  APL. 

Chapter  IV  demonstrates  the  application  of  the  complex 
LMS  algorithm  to  a  pulse  communication  problem.  The  inte¬ 
gration  time  in  this  case  is  much  shorter  than  that  of  the 
passive  sonar  case.  Two  types  of  pulse  waveforms  are  included, 
continuous  wave  (CW)  and  linear  frequency-modulated  (LFM) . 

This  simulation  program  is  implemented  in  VS  FORTRAN. 
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Chapter  V  concludes  this  thesis  by  identifying  further 
research  areas  and  other  possible  applications. 

Appendix  A  contains  the  derivation  of  the  complex  LMS 
algorithm.  Appendix  B  has  the  description  of  the  passive 
detection  program  implemented  in  APL.  Appendix  C  has  the 
description  of  the  pulse  communication  program  implemented 
in  VS  Fortran. 
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II.  THEORY  OF  SYSTEM  MODEL 


A.  OVERVIEW  OF  ARRAYS  [Refs.  5 ,6, 7 ,8] 

The  characteristics  of  the  array  elements  and  their 
arrangement  in  forming  the  array  determine  the  ultimate 
performance  of  an  adaptive  array  system  [Ref.  6].  Both  the 
linear  array  and  the  planar  array  are  examined  here. 

1.  Linear  Arrays  [Refs.  5,8] 

Consider  a  linear  array  that  has  M  equally  spaced, 
identical  point  source  elements  along  the  x-axis.  For  illus¬ 
tration,  Figure  2  shows  a  7-element  linear  array  with  uniform 
interelement  spacing  dx  and  a  plane  wave  arriving  at  the 
array  with  an  incident  angle  9  as  measured  from  the  array 
normal.  The  phasor  sum  of  all  elements  is: 


s  (t) 


M-l  .  , 

m=0 


(2.1) 


where  \p  is  the  phase  shift  between  the  nr  and  the  m+1 
element  for  m  =  0,1, 2, 3,  ...,  M-l. 


k 


d^  sin  9 
c 


sin  9 
~~c 


) 


(2.2) 


where : 


k 


wave  number. in  radians  per  meter 
A  =  wavelength  in  meters 

f  =  frequency  in  hertz 

=  interelement  spacing  in  meters 
9  *  incident  angle  measured  from  array  normal 

c  =  speed  of  wave  propagation  in  meters  per  second. 


Rewriting  equation  (2.2)  by  substituting  A  =  c/f  yields 


'P 


2*f  dX  sin  9 
c  c 


(2.3) 


The  Fourier  transform  of  equation  (2.1)  is: 


S(f)  =  F{  s  (t) }  =  /  s  (t)  e~^2lTftdt 

—  00 


(2.4) 


M-l 


=  /  l  y(t)e^m^e-=>2Trftdt 

-oo  m=0 


(2.5) 


=  l  ejmiJ;  /  y  (t)e"j27Tftdt  (2.6) 

m=0  -oo 

s ( f , 0 )  =  A ( f , 8 ) Y ( f )  (2.7) 

where  Y(f)  is  the  frequency  spectrum  of  the  incident  wave  and 
A(f,0)  is  called  the  space  factor  or  array  factor.  The  array 


factor  A(f,9)  determines  the  directional  plane  of  the  array 
in  a  plane  containing  the  array.  The  dependence  of  A(f,9)  on 
frequency,  speed  of  propagation,  element  separation,  number 
of  elements,  and  incident  angle  can  be  shown  by  rewriting 
A (f , 9 )  as : 


,2nf  ,  . 

M— 1  .  ,  M-l  3m(—r  dXSln  9) 

A(f,e)  »  l  -  l  e  c 

m=0  m=0 


(2.8) 


Summing  equation  (2.8)  yields: 


.M,  .  M  . 

sin 

A  ( f ,  9 )  =  e  ( - --)  (2.9) 

sin  j 

The  normalized  directional  pattern  is  given  by: 

2 

G  ( f ,  6  )  =  10  log  10  {  (2.10) 

M 

For  nonisotropic  (non-point  source)  elements,  it  is  necessary 
to  introduce  an  additional  factor  E(f,9)  in  equation  (2.9) 
to  include  the  directional  response  pattern  introduced  by 
each  sensor  element  [Ref.  6].  The  overall  directivity  pat¬ 
tern  then  is  given  by  the  produce  of  the  array  factor  and  the 
element  factor  [Ref.  5] .  However,  if  the  size  of  the  individual 


elements  are  small  compared  to  a  wavelength,  they  can  be 
assumed  to  be  omnidirectional  point  sources,  i.e.,  E(f,0)  =  1. 
The  effects  of  increasing  the  number  of  elements  while  main¬ 
taining  the  same  element  spacing  are  shown  in  Figures  3  and 
4  [Ref.  6] .  It  can  be  seen  that  the  main  lob  beamwidth  de¬ 
creases  as  the  number  of  elements  increases,  and  the  number 
of  sidelobes  and  nulls  increases.  In  Figures  5-8,  the  number 
of  elements  is  held  constant  at  M  =  7  while  the  interelement 
spacing  is  varied  to  illustrate  the  effects  of  elements 
spacing  on  the  directivity  pattern. 

Beamsteering  (Ref.  5]  is  accomplished  by  applying  a 
linea-r  phase  shift  across  the  line  array  as  shown  in  Figure 
9  [Ref.  6] .  The  effect  of  the  insertion  of  this  sequence  of 
phase  shifts  is  that  the  main  lobe  is  steered  to  an  angle  as 
measured  off  the  boresight  equal  to  0  where 


0 


sin^f1 


7  a. 


5} 


(2.11) 


and  o  is  the  phase  shift  between  adjacent  elements.  Figure 
10  shows  the  directivity  pattern  of  a  steered  linear  array. 

2.  Planar  Arrays  [Refs.  5,6] 

Much  of  the  analysis  done  in  linear  arrays  can  be 
extended  to  the  case  of  a  rectangular-shaped  planar  array. 

A  circular  planar  array  or  a  spherical  volume  array  would  re¬ 
quire  the  use  of  polar  and  spherical  coordinates  respectively. 
However,  array  theory  is  invariant  under  coordinate 


transformations . 
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A  planar  array  has  the  advantage  of  resolving  the 
azimuthal  and  the  elevation  angles  of  arrival  of  an  incident 
plane  wave  [Refs.  5,9].  Consider  a  rectangular-shaped 
planar  array  as  shown  in  Figure  11;  sensor  elements  are 
arranged  in  a  rectangular  grid  in  the  x-y  plane.  The  center 
of  the  array  is  usually  chosen  as  the  coordinate  origin. 

The  entire  array  has  M  elements  in  the  x-direction  with  uni¬ 
form  spacing  dx  and  N  elements  in  the  y-direction  with  uni^- 
form  spacing  dy.  The  elements  are  assumed  to  be  point  sources 
The  phasor  sum  of  the  entire  array  can  be  written  as: 

jmto  jnii; 

•'  s(t)  =  l  l  y(t)  e  X  e  y  (2.12) 

m  n 


where 


^  =  2tt(-£-)  sin9  cost); 


(2.13) 


ip  =  2tt  (— ~)  sine  sin<£ 


(2.14: 


The  directivity  pattern  of  the  planar  array  is  given  by: 


j  imp  jnijj 

A(f,6,<t>)  =  l  le  Xe  y  =  A  ( f ,  6 ,  <j> )  A  ( f ,  6 ,  <J> )  (2.15) 

x  y 


m  n 


It  follows  from  the  model  given  by  equation  (2.12)  that  the 
planar  array  beam  pattern  is  the  product  of  the  array  factors 


of  two  linear  arrays.  However,  separability  of  the  two- 
dimensional  beam  pattern  is  not  necessary  to  ensure  the 
proper  operation  of  the  LMS  adaptive  algorithm  for  a  planar 
array.  Beamsteering  is  accomplished  by  applying  appropriate 
linear  phase  shifts  for  the  row  and  column  elements.  For 
the  rectangular  array  case  under  investigation  here,  it  is 
more  convenient  to  transform  the  elevation-azimuth  (9,<j>) 
space  to  a  rectilinear  coordinate  space  (u,v)  by  the 
transformation : 

u  =  sin  0  cos  <j>  (2.16) 

v  =  sin  0  sin  <J>  (2.17) 

The  parameters  u  and  v  are  the  direction  cosines  with 
respect  to  the  x  and  y  axes,  respectively.  Figure  12  shows 
the  alternate  diagrams  for  presenting  two-dimensional  array 
beam  patterns  [Ref.  6] .  The  ranges  of  the  spherical  angles 
are  0  <  0  <  tt/2  and  0  <_  4>  <_  2tt  whereas  the  ranges  of  the 
rectilinear  coordinate  system  are  -1  <_  u  <_  1  and  -1  <_  v  <_  1 . 

B.  THE  PLANE  WAVE  MODEL 

1 .  Far-Field  Condition 

The  LMS  adaptive  filter  designed  in  this  thesis 
will  provide  spatial  resolution  for  a  planar  array  for  an 
incident  plane-wave  field.  The  plane  wave  assumption  is 
justifiable  for  a  radiating  source  located  in  the  far-field 


[Refs.  5,10] ,  due  to  wavefront  expansion.  The  far-field 
range  for  a  planar  array  is  given  by  [Ref.  5] : 


R  > 


(2.18) 


where ; 


L  0  L  0  1/2 
D  .  {(-i)2  ♦  (-£)2> 


(2.19) 


represents  the  maximum  radial  extent  of  the  transducer  array 
and  L  and  L  are  the  dimensions  of  the  planar  array  in 

x  y 

the  x  and  y  directions,  respectively. 

2 .  Propagation  of  a  Plane  Wave  from  a  Far-Field  Source 
The  plane  wave  solution  of  the  Helmholtz  wave  equa¬ 
tion  has  the  form: 


y(t,r)  =  Ae 


j  (27rft+k*r) 


(2.20) 


where  y(t,r)  is  called  the  velocity  potential,  f  is  the 
frequency,  k  is  the  propagation  vector,  and  r  the  position 
vector.  In  rectangular  coordinates, 


k  =  kx  +  ky  +  kz 

—  x  y  2 


(2.21) 


r  =  xx  +  yy  +  zz 


(2.22) 


For  a  planar  array  located  at  some  reference  location  z  =  0, 
the  velocity  potential  is: 

j  [ 2 tt £ t  +  (k  x+k  y)  ] 

y(t,r)  =  y (t ,x,y )  =  Ae  y  (2.24) 

For  a  planar  array  with  discrete  sensor  elements  located 

uniformly  in  the  x-y  plane  with  spacings  dv  and  d  respec- 

tively,  the  continuous  space  variable  x  can  be  replaced  by 

mdx  and  y  replaced  by  ndy.  If  the  time  signal  is  digitized 

for  computer  processing,  then  the  time  variable  t  can  be 

replaced  by  £T  ,  where  £  is  the  discrete  time  index  and  T 
s  s 

the  sampling  interval.  To  summarize: 

t  -  £Tg  (2.25a) 

x  -  mdx  (2.25b) 

y  -*•  ndY  (2.25c) 

The  corresponding  discrete  time,  sampled  space  signal  is 
given  by: 

j2irf£T  j  (md  k  +nd  k  ) 

y (£Ts,mdx,ndy)  =  Ae  se  x  *  (2.26) 

where  y(£T  ,md  ,nd  )  is  usually  shorted  to  y(£,m,n).  The 
propagation  vector  k  at  z  =  0  has  two  components  k  and  k 


that  can  be  related  to  the  direction  cosines  u  and  v  via 
[Ref.  5] : 


2ttu 

X 


C 


and 


2  7T  V 

A 


(2 


Substituting  equation  (2.27)  into  equation  (2.26)  yields: 


yU,m,n)  = 


j  2-rrf  £T  j~(umdx+vndY) 
Ae  se 


(2 


Let 


j2irf  £T 

y  ( £ )  =  Ae  S  ( 

represent  the  time  dependence  of  the  signal.  Equation  (2 
then  becomes : 


j~(umdY+vnd  ) 
yU,m,n)  =  y(J Ue  A  *  * 


( 


•  27a) 

.27b) 

.28) 

.29) 

28) 

.30) 


From  equation  (2.30),  it  is  easy  to  see  that  the  signal  at 
each  element  location  (m,n)  has  the  same  time  dependence 
but  has  a  different  phase  due  to  the  element  location  and 


the  direction  cosines  associated  with  the  incident  angle  of 
the  plane  wave  upon  the  array.  The  exponential  relationship 
of  the  phase  in  equation  (2.30)  also  suggests  that  the 
phases  in  the  x  and  y  directions  are  separable. 


GO 


S  (£) 


(2.32) 


j  ~ ;  umd Y+vnd  ) 

=  yU)  l  l  cd(m,n)e  '  "  1 

m  n 


If 


cd(m,n)  =  e 


2  7T 

-  j  —  ( umdx+vnd  ) 


(2.33! 


then  the  quantity  inside  the  summations  becomes  unity  and: 


sU)  =  yU)  l  l  (1) 
m  n 


(2.34) 


or 


sU)  =  MN  yU)  (2.35) 

where  s(l)  is  the  sum  over  all  elements  and  equation  (2.35) 
is  the  maximum  signal  level  possible.  This  maximum  level  is 
achieved  by  tuning  M  *N  adaptive  weights  cd(m,n)  to  conform 
to  equation  (2.33).  The  same  phase  weighting  procedure  is 
also  true  in  the  frequency  domain,  in  fact,  the  implementation 
of  phase  weighting  is  inherently  a  frequency  domain  opera¬ 
tion.  Since  phase  weight  equation  (2.33)  is  a  function  of 
wavelength  A  and  A  =  c/f,  the  proper  phase  weight  for  co¬ 
phasing  at  each  element  is  a  function  of  frequency  f.  Thus 
for  each  valid  frequency  component  in  the  signal,  a  different 


36 


set  of  phase  weights  {cd(m,n)}  must  be  generated.  Consider 
the  DFT  with  respect  to  time  of  equation  (2.32): 

j^(umdy+vnd  ) 

S(q)  =  Y (q)  l  l  cd (m,n) e  x  (2.36) 

m  n 


\ 

i 


where  q  is  the  frequency  index  and  Y(q) 
If  equation  (2.33)  holds,  then: 


is  the  DFT  of  y (1). 


S  (q)  =  Y (q)  I  I (1)  =  MN  Y(q) 

m  n 


(2.37) 


Only  valid  spectral  lines  will  be  processed. 

2 .  The  Frequency  Domain  LMS  Adaptive  Filter 

The  general  frequency  domain  LMS  adaptive  algorithm 
is  derived  in  Appendix  A.  Suppose  that  the  time  signal  zU) 

A 

is  the  reference  or  desired  signal  and  z^(£)  is  the  normalized 
sum  of  all  signals  in  the  planar  array.  In  the  frequency 
domain,  the  reference  signal  in  the  q  DFT  bin  is: 

L-l  -j^r£q 

Z(q)  =  l  z(£)e  L  (2.38) 

£=0 


and  the  estimated  output  in  the  frequency  domain  is: 


Substituting  equation  (2.31)  with  U)  =■  yields : 

L-l  -j^£q 

Zi(q)  =  l  mn  Z  I  cd.(m,n)y(£,m,n)e  (2.40) 

£=0  m  n 

where : 

£  is  the  time  index;  £  =  0,1,..., L-l 

q  is  the  DFT  bin  index;  q  =  0,1,..., Q-l 

i  is  the  complex  phase  weight  iteration  number. 

The  DFT  operation  with  respect  to  time  in  equation  (2.40) 
can  be  performed  first  to  yield: 

Zi(q)  =  I  1  cdi(m,n)Y(q,m,n)  (2.41) 

m  n 

where : 

L-l  “ 

Y (q,m,n)  =  [  y(£,m,n)e  (2.42) 

£  =  0 

Based  on  the  complex  LMS  algorithm  [Ref.  1] ,  the  adaptive 
filter  output  in  the  q  bin  is  given  by  equation  (2.41). 

The  error  signal  is  generated  by  comparing  the  desired 
(reference)  signal  to  the  adaptive  filter  output  (estimate) . 
The  error  is  denoted  by  e.,  where: 


The  estimate,  equation  (2.41),  is  formed  by  applying  the 

till 

phase  weights  (cd^(m,n)}  of  the  i  iteration  to  each  element 
in  the  planar  array.  The  complex  weight  cd^(m,n)  is  updated 
recursively  as  follows: 


cdi+i  (m,n)  =  cd^n^n)  +  2uie1Y*  (q,m,n) 


(2.44) 


where : 

m  =  0,1,..., M-l 

n  =  0,1,... , N-l 

(*)  denotes  complex  conjugate 

U ■  =  feedback  coefficient,  a  parameter  that  con- 

1  trols  the  rate  of  convergence,  algorithm 

noise,  and  the  stability  of  the  algorithm 
[Ref.  4]  . 

til 

From  equation  (2.44),  it  can  be  seen  that  the  i+1  weight 
cdi+^(m,n)  may  have  magnitudes  larger  than  unity.  This 
growth  in  magnitude  is  undesirable  for  the  purpose  of 
spatial  resolution  since  spatial  resolution  depends  on  the 
relative  phase  between  adjacent  elements  to  resolve  the 
direction  of  wave  arrival.  Thus,  a  normalization  is  neces¬ 
sary  to  bring  equation  (2.44)  back  to  unity.  This  normali¬ 


zation  is: 


(2.45) 


cdi+l  (m,n)  cdi+1(m,n)/|cdi+1(m,n)  | 


These  updated  and  normalized  phase  weights  can  now  be  applied 
to  equations  (2.41),  (2.43),  and  (2.44)  in  sequence  to 
compute  the  next  set  of  phase  weights.  This  iterative 
process  stops  when  predetermined  criteria  are  met.  At  that 
point,  the  set  of  phase  weights  {cd^(m,n)}  can  be  used  to 
find  the  direction  cosines  of  the  incident  plane  wave. 
However,  the  phase  angles  of  the  phase  weights  {cd^(m,n)} 
may  have  been  wrapped  around  an  integer  multiple  of  2tt.  The 
procedure  for  phase  unwrapping  is  explained  in  Sections  4  and 
5  of  this  chapter. 

The  feedback  coefficient  in  equation  (2.44)  con¬ 
trols  the  rate  of  convergence  and  the  stability  of  the  LMS 
adaptive  filter.  Robbins  and  Monro  [Ref.  11]  showed  that 
the  adaptive  weights  (cd^(m,n)}  will  converge  to  the  optimum 
result  if  yu  is  allowed  to  decrease  with  the  iteration  index 


< 


00 


00 


l 

i=l 


2 

Ui 


A  coefficient  y^  that  satisfies  the  above  conditions  will 
work  as  long  as  the  signal  and  noise  inputs  are  truly 
stationary  but  will  not  be  satisfactory  for  a  filter  operat¬ 
ing  in  a  slowly  varying  environment.  Widrow's  LMS  algorithm 
[Ref.  1]  uses  a  constant  value  of  y  satisfying  the  inequality 


where  X  is  the  largest  eigenvalue  of  the  correlation 
max 

matrix  of  the  input.  Although  this  matrix  is  typically  not 
known  a  priori,  some  bound  can  be  set  up  by  examining  equation 
(2.44).  If  stationarity  can  be  assumed,  it  is  possible  to 
update  y^  every  iteration  in  order  to  obtain  the  optimum  set 
of  phase  weights.  In  Appendix  A  a  simple  method  of  updating 
the  feedback  coefficient  y^  is  proposed  to  improve  the  per¬ 
formance  of  the  LMS  adaptive  filter. 

3 .  Applying  the  Frequency  Domain  LMS  Adaptive  Filter 
to  a  Planar  Array 

Given  a  rectangular  planar  array  with  M  *N  elements, 
there  are  several  ways  to  process  the  signal  from  the  array 
to  achieve  spatial  resolution  in  both  azimuth  and  elevation. 
Three  different  ways  are  considered  in  this  thesis. 
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a.  Orthogonal  Linear  Arrays 

Two-dimensional  spatial  resolution  is  possible 
by  just  considering  one  linear  array  in  the  x-direction  and 
one  linear  array  in  the  y-direction.  A  total  of  M+N-l  ele¬ 
ments  out  of  M  x n  are  used.  This  scheme  is  useful  when 
processing  time  is  limited.  Figure  14  illustrates  the  choice 
of  the  center  linear  arrays  for  this  scheme.  However,  any 
two  orthogonal  linear  arrays  in  the  planar  array  can  be 
used.  The  recursive  equations  needed  to  implement  this 
algorithm  are  divided  into  two  sets;  one  set  for  the  linear 
array  in  the  x-direction  and  the  other  for  the  y-direction. 

In  the  x-direction  the  estimate  is: 

.  .  M— 1 

Zx  (q)  =  M  ^  Vnt)Y(q,m,n  -nQ)  '  '  (2.46) 

i  m=0 

where : 

n  =  constant  y-direction  index 
o 

c^(m)  =  unity  magnitude  phase  weight. 

The  error  is  the  difference  between  the  reference  and  the 
estimate. 

e  =  Z (q )  -  Zv  (q)  (2.47) 

1  1 

The  recursive  update  for  phase  weights  is: 


ci+l  (m)  =  Ci(m)  +  2ux  ex  Y*(q,m,n=nQ) 

i  i 


(2.48) 


The  phase  of  the  new  update  is : 


ci+l  (in)  *■  ci+1(m)/|ci+1(m) 


In  the  y-direction,  the  procedure  is  similar; 


i  N-l 

Estimate:  Z  (q)  =  -  [  d .  (n)  Y  (a  ,m  =m  ,n) 

yi  N  n=0  1  '  ° 


Error:  e  =  Z(z)  -  Z  (q) 

yi  yi 


(2.49) 


(2.50) 


(2.51) 


Update:  di+l  35  (n)  +  2u  e  Y*(q,m=m  ,n)  (2.52) 

yi  yi  ° 


Normalization:  d^+^ (n) 


di+i(n)/|di+l(n) 


(2.53! 


The  convergence  constants  ux  and  are  usually  set  to  be 
equal  since  the  statistics  in  the  orthogonal  directions  of 
a  planar  array  can  be  assumed  to  be  the  same  for  the  obser~ 
vation  time  of  most  systems. 

b.  Two-dimensional  Array 

This  scheme  uses  all  M  *N  elements  and  therefore 
it  can  realize  the  full  array  gain  of  equation  (2.35) .  The 


phase  weights  cd(m,n)  are  not  assumed  to  be  separable.  The 
equations  for  the  LMS  adaptive  filter  in  the  frequency 
domain  are: 


Estimate:  (q)  =  5^  [  [  cd^ (m,n) Y (q,m,n)  (2.54) 

m  n 

Error:  e±  =  2 (q)  -  Z± (q)  (2.55) * 

Update:  cd^+^(m,n)  =  cd^(m,n)  +  2u^e^Y* (q,m,n)  (2.56) 

Normalization:  cdi+1(m,n)  cd^+1  (m,n)  /  |  cdi+1  (m,n)  |  (2.57) 


c.  Separable  Two-dimensional  Array 

As  mentioned  in  the  discussion  on  planar  arrays 
and  plane  waves,  the  form  of  a  plane  wave  suggests  that  the 
phase  of  a  signal  at  an  element  (m,n)  is  separable.  This 
scheme  then  uses  the  separability  property 

cd(m,n)  =  c(m)d(n)  (2.58) 

to  implement  a  two-dimensional  LMS  adaptive  filter.  All 
MxN  elements  are  used  but  only  M+N  phase  weights  need  to  be 
updated  recursively.  The  equations  are: 

Estimate:  Z^(q)  =  ^  £  d.(n)  l  c.  (m) Y (q,m,n)  (2.59) 

n  1  m  1 
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A 

Error:  =  Z(q)  -  Z^(q)  (2.60) 

Updates:  ci+i  =  +  2uiej.(I  d^ (n) Y (q,m,n) )  *  (2.61) 

n 

di+i (n)  =  di(n)  +  2uiei([  (m) Y (q,m,n) ) *  (2.62) 

“  m 

Normalizations:  ci+]_(m)  *■  ci+i  /  i  ci+i  (m>  !  •'  (2.63) 

di+l(n)  "  di+i^)/ldi+i(nH  (2.64) 

All  three  of  the  aforementioned  schmes  are  implemented  and 
their  results  compared.  At  the  start  of  all  three  algorithms, 
the  initial  phase  weights  are  set  to  the  boresight  of  the 
planar  array,  i.e.,  magnitude  equal  to  unity  and  phase  angle 
to  zero.  The  normalization  of  phase  weights  to  unity  forces 
the  spatial  transfer  function  of  the  LMS  adaptive  filter  to 
have  unit  magnitude.  The  steady  state  phase  response  is 
designed  to  phase  align  all  sensor  elements  in  an  element- 
by-element  fashion.  More  discussion  on  this  topic  can  be 
found  in  Appendix  A. 

A 

4 .  Extracting  Estimates  of  the  Direction  Cosines  u  and 
v  from  Phase  Weights 

A  A 

To  extract  u  and  v  from  the  orthogonal  linear  arrays 

and  the  separable  two-dimensional  array  cases  discussed  in  ; 

Section  C.3.a.  and  C.3.a.c.,  only  M  elements  in  the  x-direc- 

tion  and  N  elements  in  the  y-direction  need  to  be  considered. 

' 

I 
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a.  Direction  Cosine  Estimates  for  Linear  Arrays 

Consider  that  in  cases  3a.  and  3c.,  the  phase 

weights  c (m)  ,  d(n)  have  reached  a  steady  state.  The 

objective  at  this  point  is  to  relate  the  phase  angles  of 

these  two  sets  of  phase  weights  to  their  respective  direction 

cosines.  Let  £  (m)  be  the  phase  angle  of  c (m)  and  £  (n)  be 
x  y 

the  phase  angle  of  d(n),  i.e., 


,  _  ^__-l  r  lm[c  (m)  ]  , 

£x(m)  -  tan 


(2.65a) 


/_v  _  ^__-l  f  Im[d  (n)  ]  ! 

Vn)  "  ReT'd  (n)  ] 


(2.65b) 


It  can  be  seen  from  equation  (2.33)  and  using  the  concept 
of  separability  that: 


5x<m>  * 


2tt  . "  .  , 

-  — (u  mdx) 


(2.66a) 


£y(n)  = 


2tt  ,  ~  .  , 

'  IT (V  n  dY 


(2.66b) 


Solving  equation  (2.66)  for  u  and  v  yields 


~ A  E  (m) 


u  = 


,  m  ?  0 


(2.67a) 


“A£  (n) 


,  n  ^  0 


(2.67b) 


MNiunuiuiiuniffUffiM  mail  rnmninmiTfi  inn  iiipuuwui»i«  «* 


Knwj'.  ^  v  ktjuv 


*mw.iuwNPwnn  u  ivmiimpiwwi 


Thus,  in  the  x-direction,  there  are  M-l  estimates  of  u; 
while'  in  the  y-direction,  there  are  N-l  estimates  of  v. 

To  find  an  estimate  of  the  direction  cosines  from  equation 
(2.67),  one  needs  to  take  an  arithmetic  average  of  possible 
estimates : 


u  - 


M-l 


M-l  -X£  (m) 

l 


m=l 


2iTmd. 


(2.68a) 


v  = 


!  N;1  (n) 


N-. 


I 

n=l 


TimcT 


(2.68b) 


Equations  (2.68a)  and  (2.68b)  will  be  referred  to  as  the 

/V  /N 

point  by  point  method.  Another  way  of  finding  u  and  v  makes 
use  of  linear  regression  [Ref.  13],  Consider  equation  (2.66) 


where  £x(m)  is  a  linear  function  of  m  with  slope  equal  to 


2  7T  ^ 

— ;j-udx  and  C  (n)  is  a  linear  function  of  n  with  slope  equal 


2tT 


to  — — vd^.  Usin<?  a  linear  regression  fit  of  M  data  points 
vs.  the  element  number  m,  the  slope  and  intercept  of  the  line 


£x(m)  can  be  calculated.  The  same  procedure  can  be  used 


for  £y(n).  Let  the  slope  in  the  x-direction  be  s^  and  the 


slope  in  the  y-direction  be  S2*  Then: 


2tt 

-  T  udx 


(2.69a) 


2  it  " 

"  "T  vdY 


(2.69b) 
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Thus,  u  and  v  can  be  solved  by  rearranging  equation  (2.69) 
to  the  form: 


u  = 


(2.70a) 


V  '  2^ 


(2.70b) 


The  estimates  u  and  v  obtained  from  the  linear  regression 
method  represent  the  best  linear  least-squares  fit  of  the 
observed  data. 

b.  Direction  Cosine  Estimates  for  Planar  Arrays 
The  two-dimensional  array  discussed  in  Section 
3.b  of  this  chapter  does  not  require  the  phase  weights  to 
be  separable.  Consider  that  the  set  of  phase  weights 
{cd(m,n)}  has  reached  a  steady  state.  Let  £  (m,n)  be  the 

xy 

phase  angle  associated  with  the  cd(m,n).  Then: 


5xy(m,n) 


.  (m,n)  ]  ■, 

tan  Wca^HTTTjr 


(2.71) 


From  equation  (2.33),  it  can  be  seen  that 


£  (m,n)  = 

xy 


-  — ^-(u  mdx  +  v  ndy) 


(2.72) 


In  general,  £x^(m,n)  can  be  a  more  complicated  function  of  m 
and  n.  For  instance,  in  the  near-field  problem,  one  needs  to 


WWW 


modify  equation  (2.72)  to  contain  quadratic  phase  terms  to 
account  for  wavefront  curvature  [Ref.  5].  For  the  plane 
wave  model,  equation  (2.72)  adequately  describes  the  phase 
weights  needed  to  steer  the  directivity  pattern  of  the  planar 

A  A 

array  to  the  direction  corresponding  to  u  and  v. 

The  point  by  point  method  is  applicable  here 

A  /A 

to  find  u  and  v  given  a  steady  state  phase  angle.  However, 
a  linear  regression  fit  of  £Xy(m,n)  vs.  the  element  index 
m  and  n  appears  to  be  more  suited  to  this  problem.  Rewriting 
equation  (2.72)  in  the  form  of  the  equation  for  a  plane  in 
three-dimensional  space  yields: 

£xy(m,n)  =  (I^lGdx)m  +  (2.73) 

Equation  (2.73)  describes  a  plane  with  slope  -yudx  in  the 

2tt  a 

m-direction  (x-direction  and  slope  — —  vdy  in  the  n-direction 
(y-direction) .  Again,  let  the  slopes  be: 

S1  =  "  X  ^  dX  (2.74a) 

s2  =  -  ~  v  dy  (2.74b) 

Thus,  equation  (2.74)  is  identical  to  equation  (2.70)  and 

A  /\ 

the  direction  cosine  estimates  u  and  v  are  found  by  equation 
(2.71) ,  i.e. , 
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I 


(2.71a) 


(2.71b) 


This  result  is  not  surprising  since  the  exponential  represen¬ 


tation  of  the  plane  wave,  equation  (2.20),  is  inherently 


separable . 


5.  Unwrapping  the  Steering  Phase  Weights 


a.  Linear  Array  Unwrapping 


The  proper  phase  weights  for  beamsteering  are 


given  by  equation  (2.66): 


?x(m)  = 


-2it  ,  ~  .  . 

— (u  mdx> 


(2.66a) 


V": 


-2tT~  ,  . 

(v  ndy) 


(2.66b) 


Consider  a  7-element  linear  array  lying  in  the  x-direction 


with  the  center  element  as  the  reference  element.  The  element 


index  then  runs  from  m  =  -3,  -2,  -1,  0,  1,  2,  3.  If  u  is 


equal  to  0.55  and  dv  =  A/2,  then  equation  (2.67a)  reduces  to: 


(m)  =  -  Hum 


(2.75) 


=  -0 . 55  tt  m 


■* 


The  following  table  shows  the  required  phase  weights  E  (m) 

A 

A 

needed  to  steer  the  beam  to  u  =  0.55. 


TABLE  1 

PHASE  WEIGHTS  FOR  BEAMSTEERING 


m 

-3 

-2 

-1 

0 

1 

2 

3 

1.65w 

l.lTT 

.55* 

0 

-.5577 

-1.177 

-1.6  577 

q(m) 

-.3571- 

-.977 

.5577 

0 

-.5577 

0.9tt 

.3577 

The  phase 

factors 

j?x(m) 

le 

}  are  a 

set  of 

complex 

numbers 

in 

the  complex  plane.  The  angle  of  any  complex  number  must  lie 
within  a  2tt  interval.  The  interval  chosen  here  is  [-17,77]. 
This  means  that  any  angle  £  (m)  that  is  outside  the  range 

A 

[ — -rr , tt ]  will  be  wrapped  around  to  an  angle  £^(m)  that  is 
inside  the  range.  This  property  can  be  shown  as  follows: 


ej(S+2irk)  =  e j  6e j 27Tk  =  0j3 


since 

ei2irk  =  i  for  k  ■  0,±1,  ±2, . . . 


therefore 


4  O 

eJ  has  modulo  2 tt. 
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Referring  back  to  Table  1,  the  observed  angles  £'(m)  are 

X 

wrapped.  If  no  processing  is  done  to  unwrap  the  observed 
angles,  the  equations  derived  in  Section  C.4  of  this  chapter 
will  not  apply  for  all  permissible  values  of  u  and  v.  For 

r-.  /V 

small  values  of  u  or  v,  the  needed  phase  weights  do  not  wrap 
around  but  the  spatial  range  of  interest  is  severely  restricted 
[Ref.  4].  In  tracking  systems,  the  above  restriction  in  look 
direction  can  be  justified  since  a  crude  target  direction  is 
usually  provided  by  a  search  array.  The  maximum  spatial 
window  of  an  M  element  linear  array  without  the  unwrapping 
of  the  steering  phase  weights  can  be  calculated.  For  example, 
given-  that  the  element  spacing  is  dx  =  X/2,  and  the  maximum 
permissible  magnitude  of  E,  (m)  is  tt  in  the  range  [-it, it],  then 
equation  (2.66a)  becomes: 


£  (m)  |  =  tt  = 

sx  'max 


(2.77) 


or 


In  general. 


A 


(2.78) 


where  M  is  the  number  of  elements  in  the  array.  The  corres¬ 
ponding  angular  coverage  in  terms  of  the  polar  coordinate 
angle  9  is  (for  0=0): 

0  =  sin’1(l“!max>  (2-79) 

I 

| 

;  The  total  angular  coverage  is  29.  For  the  case  of  the  7- 

I  element  linear  array,  this  corresponds  to  a  coverage  of  about 

I 

40°  out  of  a  range  of  180°.  Figure  15  shows  the  expected 
phase  angles  required  for  beamsteering  vs.  element  number. 

|  Figure  16  shows  the  wrapped  angles  vs.  element  number.  It 

;  should  be  noted  that  either  set  of  these  angles  (phase 

I 

|  weights)  will  steer  the  beam  to  the  proper  direction.  The 

difficulty  with  the  wrapped  (observed)  angles  is  that  the 

A  /V 

direction  cosines  u  and  v  cannot  be  directly  estimated  using 
the  methods  developed  in  Section  C.4  of  this  chapter.  In 
order  to  unwrap  the  observed  angles  in  Table  1,  consider 
equation  (2.76).  It  can  be  seen  that  the  observed  angles 
differ  from  the  angles  generated  from  equation  (2.66)  by 
an  amount  of  ±2-rrk  where  k  =  0,1,2,....  In  order  to  unwrap 
the  observed  phase  angles,  it  is  necessary  to  find  out  which 
elements'  phase  angles  have  been  wrapped  around  and  by  what 


RE 


ASE 


integer  multiple  of  2*.  Let  w(m)  be  the  unwrap  factor  such 
that : 


5x(m)  =  Sx(m)  =  £^(m)  +  w(m)  (2.80) 

A 

Table  2  shows  (for  u  =  0.55  and  =  A/2)  how  the  implemen¬ 
tation  of  equation  (2.80)  will  yield  the  phase  weights 
required  for  spatial  resolution. 


TABLE  2 

PHASE  UNWRAPPING 


m 

-3 

-2 

1 

0 

1 

2 

3 

Cx(m) 

1.65rr 

1.1* 

.55* 

0* 

-.55* 

-1.1* 

-1.65* 

^<m> 

-.35* 

-.9* 

.55* 

0 

-.55* 

0.9* 

0.35* 

w(m) 

+2* 

+2* 

0 

0 

0 

-2* 

-2* 

S>) 

1.65* 

1.1* 

.55* 

0 

-.55* 

-1.1* 

-1.65* 

Observe  from 

Table  2  and  equation  (2. 

66)  that  the  phase 

angles  for  elements  m  = 

=  -1  and 

m  =  1  i 

do  not  wrap  around  as 

long  as  dx  £ 

A 

2  ' 

Recall  also  that  the 

center 

element 

is 

chosen  as  the  reference  element.  It  is  possible  then  after 
a  set  of  steady  phase  weights  has  been  computed  that  an 
estimate  of  the  direction  cosine  can  be  computed  using  the 
phase  angles  at  elements  m  =  -1  and  m  =  1  in  equation  (2.67) 


below. 


u  <m)  = 


-X5x(m) 


(2.67) 


-  A/2, 


u  (m) 


-5x(m) 


(2.81) 


The  estimate  of  direction  cosine  u  using  only  the  information 
(phase  angles)  at  m  =  ±1  is  denoted  as  where 


u  =  -rtu  (1)  +  u  (-1)  ] 

g  2 


(2.82) 


Equation  (2.82)  will  yield  a  good  estimate  of  u  as  long  as 
dX  —  J'  The  next  steP  i-n  this  process  of  phase  unwrapping 
is  to  use  the  result  from  equation  (2.82)  to  generate  a  set 
of  projected  phase  weights  {£x(m)}  using  the  following 
relation: 


«*'»>  - 


-2TTmdvu 

A  C 


(2.83) 


Recall  from  equation  (2.76)  that  a  phase  angle  outside  the 
rangle  is  mapped  to  an  angle  within  {-it, it].  By 

examining  the  magnitude  of  the  projected  angle,  it  is  possible 
to  determine  how  many  multiples  of  2tt  were  lost  due  to  the 

A 

modulo  2 tt  property  of  complex  numbers.  The  sign  of  ug  and 


the  location  of  the  element  m  determine  the  sign  of  the 
unwrap  factor  w(m) .  For  linear  arrays  with  an  odd  number  of 
elements,  the  unwrap  procedure  can  be  described  as  follows: 


A 

if  ( k—  2 )  it  <_  |£(m)|  <  kir 


(2.84) 


and  if  | £ (m) |  does  not  lie  between  the  limits  in  equation 
(2.84),  the  value  of  k  is  decremented  by  two  and  the  inequality, 

A 

equation  (2.84),  is  tested  again  until  1 ^ (m) |  falls  within  a 
2it  interval,  then 


| w(m)  |  =  (k-l)ir  (2.85) 

For  an  M-element  (odd)  array,  the  initial  value  of  k  is 
M-l 

(— ) .  The  sign  of  w(m)  is  determined  by: 

A 

sign  w(m)  =  -  sign  (m)  sign  (u^)  (2.86) 

where : 

m  n  /M-l, 

m  (  ^  //••.,  2,  l , 0 , 1 ,  .  .  .  ,  (  2  )  i 

i.e.,  using  the  center  element  as  the  reference.  A  similar 
procedure  can  be  utilized  in  an  array  with  an  even  number  of 
elements.  The  reference  used  should  then  be  the  point  be¬ 
tween  the  two  center  elements  since  there  is  no  need  for  the 

t 

I 

I 

) 


reference  point  to  coincide  with  the  location  of  a  sensor 
element.  This  choice  of  reference  takes  full  advantage  of 
the  resulting  symmetry.  From  Table  2,  it  can  be  seen  that 
only  half  of.  the  elements  need  to  be  examined  using  equations 
(2 . 84) - (2 . 86) .  The  unwrap  factors  {w(m)}  for  the  other  half 
are  simply  the  negative  of  the  first  half.  This  unwrap 
procedure  is  good  for  the  two  array  configurations  discussed 
in  Sections  3a.  and  3c.  of  this  chapter.  The  unwrap  proce¬ 
dure  for  the  two-dimensional  array  is  similar  but  the  compu¬ 
tation  is  a  little  more  involved.  The  unwrapping  procedure 
for  the  linear  array  in  the  y-direction  is  identical  with 

the  substitution  d  -*•  d  ,  m  -*•  n,  £  (m)  -*•  £  (n),...,  etc. 

AY  x  y 

b.  Two-dimensional  Array  Unwrapping 

The  proper  phase  weights  to  steer  a  beam  to 
(u,v)  are  given  by: 

£Xy  (m,n)  =  — jp(u  mdx  +  v  ndy)  (2.87) 

The  unwrapping  procedure  is  best  illustrated  by  considering 

A  /\ 

the  following  example.  Consider  the  case  (u,v)  =  (-0.7, +0.7) 
and  dx  =  dy  =  A/2.  Equation  (2.87)  then  becomes: 

A  A 

£  (m,n)  =  -7r(um  +  vn)  (2.88) 

xy 

=  -it  (-0 . 7m  +  0 . 7n) 


Assume  also  that  we  are  given  a  planar  array  with  M  elements 
in  the  x-direction  and  N  elements  in  the  y-direction  (the 
corresponding  element  indices  are  (x,y)  -*■  (m,n)  ,  with  both 
M  and  N  equal  to  5.  For  the  sake  of  symmetry,  let  the 
element  indices  run  as  follows: 


m 


-  (^)  ,-(^)+l . 0,1 


,N-1.  ,N-1 .  -  .  .  ,N-1. 

n  *■(  2  )  ,  ~  (  2  )"^1,*«»,0, !,••#,(  2  ) 


In  this  case,  for  M  =  5, 

m  =  -2 , -1 ,0,1,2. 


Similarly,  for  the  orthogonal  direction, 

n  =  -2, -1,0, 1,2. 

The  obvious  choice  of  reference  element  is  (m,n)  =  (0,0), 
i.e.,  the  center  element.  The  desired  phase  weights  for 
this  example  are  given  in  Table  3. 

The  phase  weights  marked  within  the  two  triangles 
in  Table  3  will  be  wrapped  around,  so  the  actual  observed 
phases  when  the  phase  weights  reach  steady  state  are  tabu¬ 
lated  in  Table  4. 
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TABLE  3 

DESIRED  PHASE  WEIGHTS  UXy(m,n)}  GIVEN 
THAT  u  =  -0.7,v  =  +0.7 


-2 . 8tt 

-2. Itt 

-1 . 4tt 

-2. Itt 

-1.4tt 

-0 . 7tt 

-1 .  4tt 

-0 .  7tt 

0 

1= 

• 

o 

1 

0 

0 . 7tt 

0 

0 . 7tt 

1 . 4it 

-2 

-1 

0 

■0 . 7tt 
0 

0 . 7it 

1 . 4  7T 

2 .  Itt 


TABLE  4 

OBSERVED  PHASES  OF  THE  STEADY  STATE 
ADAPTIVE  WEIGHTS  ^y(m,n) 


i 

2 

-0 . 8tt 

-O.lTT 

0 . 6tt 

-0 . 7tt 

0 

1 

-O.lir 

0 . 6tr 

-0 . 7tt 

0 

0 . 7tt 

t 

.« 

* 

0 

0 . 6tr 

i 

o 

• 

-4 

0 

0 . 7tt 

i 

o 

• 

<7\ 

=1 

- 

*i 

i. 

-1 

-0 . 7tt 

0 

0 . 7  it 

-0 . 6tt 

O.lTT 

'i 

»i 

»i 

-2 

0 

0 . 7tt 

-0 . 6tt 

O.lTT 

t= 

00 

• 

o 

1 

* 

• 

-2 

-1 

0 

1 

2 

Let  w(m,n)  be  the  two-dimensional  unwrap  factor,  i.e.,  the 
unwrapped  phase  is : 

£  (m,n)  =  £"  (m,n)  =  £'  (m,n)  +  w(m,n)  (2.89) 

xy  xy 

By  comparing  (5x^(m,n)}  in  Table  3  with  {£^  (m,n)}  in  Table 
4,  it  can  be  seen  that  w(m,n)  must  be  equal  to  the  tabu¬ 
lated  values  in  Table  5. 

TABLE  5 

UNWRAP  FACTORS  w(m,n) 


n 


2 

-2  TT 

-2tt 

~2tt 

0 

0 

1 

-2n 

-2  IT 

0 

0 

0 

0 

-2tt 

0 

0 

0 

2ir 

-1 

0 

0 

0 

2tt 

2tt 

-2 

0 

0 

2  TT 

2tt 

2tt 

-2 

-1 

0 

1 

2  m 

To 

generate 

w(m,n)  ,  estimates  of  the 

direction 

cosines 

A 

U 

and  v  must 

be 

computed 

using  elements 

at  m  =  ± 

1  and 

n 

=  ±1.  Let 

A 

U 

g 

A 

and  v  be 

g 

those  estimates  obtained  using 

point  by  point  method  below: 


j[u(l)  +  u(-l)] 


(2.90) 


[v  (1)  +  V (-1)  ] 


(2.91; 


For  dx  =  =  A/2,  equation  (2.88)  can  be  applied  to  com- 

A 

pute  the  projected  phases  {?x^(m,n)}  by  the  relation: 


£  (m,n)  =  — rr{u  m  +  v  n}  for  all  m  and  n  (2.92; 

xy  g  g 


The  set  of  projected  phases  (angles)  are  then  examined  to 
decide  the  proper  unwrap  factor  for  a  particule  element 
(m,n) .  The  logic  is  as  follows:  for  each  n,  all  elements  m 
are  examined. 


Check  (k-2)7T  <_  ]£Xy(m,n)j  <  kn 


(2.93) 


and  if  |£Xy(m,n)  |  does  not  lie  between  the  limits,  then  k  is 
decremented  by  two  and  equation  (2.93)  is  tested  again.  If 

A 

|£Xy(m,n) |  does  lie  between  the  limits,  then 


I  w  (m,n)  I  =  ( k— 1 )  tt 

xy 


(2.94) 


This  procedure  applies  when  M  is  odd  and  the  initial  value 
M-l 

of  k  is  (— ) •  The  sign  of  w(m,n)  is  determined  by: 


sign (w (m,n) ]  =  -sign [sign (m) sing (Ug) +sign (n) sign (v  ) ]  (2.95) 
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where : 


M— 1  M-l 

m  =  -(^±),...,-l,0,l,...,(^i) 


n  _  _  (N-l .  -I  n  ,  ,N-1. 

n  2  /  /  \  2  ' 


For  the  next  value  of  n,  equations  (2.93),  (2.94),  and  (2.95) 
are  repeated  for  all  values  of  m.  This  continues  until  the 
last  value  of  n  is  reached.  It  is  possible  though  to  examine 
only  half  of  the  planar  array  since  the  unwrap  factors  for 
the  other  half  are  the  negative  of  that  of  the  first  half. 

To  ensure  symmetry  about  the  reference  element,  it  is  neces¬ 
sary  to  rotate  the  phase  angle  at  the  reference  element  to 
zero.  This  can  be  accomplished  by  multiplying  the  phase 
weights  of  all  M  N  elements  by  the  complex  conjugate  of 
the  reference  phase  weight,  i.e.. 


cd.  (m,n)  •<-  cd.  (m,n)  cd.  (m  .n  ) 

1  1  i  o  o 


(2.96] 


where  m  and  n  are  the  indices  locating  the  reference  element. 
This  operation  will  ensure  that  the  unwrapping  procedure 
will  work  properly.  For  linear  array  phase  weights,  this 
phase  centering  should  also  be  completed  before  unwrapping. 
The  equations  are: 


ci  (m) 


c^m)  c?(mQ) 


(2.97) 
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cL  (n) 


d.(n)  d“(no) 


(2.98) 


where  mQ  is  the  index  locating  the  reference  element  of  the 


linear  array  along  the  x-direction  and  nQ  is  the  index 


locating  the  reference  element  of  the  linear  array  along  the 


y-direction. 


6 .  Summary 


The  original  complex  LMS  adaptive  filter  [Ref.  1] 


requires  three  necessary  modifications  to  make  it  useful  for 


estimating  the  direction  cosines  of  an  incident  plane  wave. 


Without  the  following  modifications,  accurate  spatial  resolu¬ 


tion  is  not  possible: 


-  normalization  of  the  adaptive  complex  weights  (phasors) 
to  unity  magnitude  after  each  iteration. 


unwrapping  of  the  observed  steady  state  phase  angles 
to  extend  the  spatial  coverage  to  the  full  range  of 
u  and  v. 


allowing  the  feedback  coefficient  y  to  decrease  with 
increasing  iteration  number  to  achieve  convergence  to 
an  optimal  set  of  phase  weights  and  to  realize  a 
robust  filter. 


D.  NOISE  MODEL 


In  the  SONAR  environment,  the  ambient  noise  field  is  a 


composite  of  many  different  noise  sources.  Therefore,  using 


the  Central  Limit  Theorem  [Ref.  13],  the  ambient  noise  can 


be  modeled  quite  adequately  as  additive  white  Gaussian  noise. 


Intentional  jamming  is  not  considered  in  this  thesis.  How¬ 


ever,  the  use  of  an  adaptive  filter  to  place  a  null  at  the 


spatial  location  of  a  jamming  signal  has  been  studied  exten¬ 
sively  [Refs.  2,4,15].  The  noise  corrupted  received  signal 
is  given  by: 

r(£,m,n)  =  y(i.,m,n)  +  n(£,m,n)  (2.99) 


where  y(£,m,n)  is  the  sampled  signal  represented  by  equation 

> 

(2.26)  and  n(Z, m,n)  is  white  Gaussian  noise  time  samples 

2 

with  zero  mean  and  noise  power  aN-  If  A  is  the  signal 
amplitude  (see  equation  (2.26)),  then  the  signal-to-noise 
ratio  is  given  by: 


or 


SNR 


(2.100) 


(SNR)dB  =  10  log 1Q (SNR)  dB  (2.101) 

The  performance  of  the  frequency  domain  LMS  adaptive  filter 
was  tested  for  various  signal-to-noise  ratios.  The  noise 
environment  is  assumed  to  be  stationary  during  the  look 
interval  of  the  SONAR  system. 


III.  LOW  FREQUENCY  PASSIVE  SONAR  TARGET  LOCALIZATION 


i 

i 


The  objectives  in  passive  SONAR  are  to  detect  and  possibly 
classify  noise  sources  in  the  ocean.  Most  of  the  noise  gener¬ 
ated  by  vessels  is  concentrated  in  the  low  frequency  range 
(30-1K  Hz) .  Acoustic  energy  in  this  frequency  range  is 
capable  of  long  range  propagation  and,  as  a  result,  most 
long  range  detection  systems  in  the  underwater  environment 
operate  at  these  low  frequencies.  Propeller  cavitation, 
machinery  noise,  and  wake  are  the  major  sources  of  such  noise. 
Under  certain  situations,  very  long  range  detection  has  been 
demonstrated  in  this  frequency  range.  In  a  long  range 
detection  scneario,  knowledge  of  the  elevation/depression 
angle  is  necessary  to  resolve  potential  range  ambiguities 
in  convergence  zone  (CZ)  problems.  Therefore,  the  use  of  a 
planar  array  is  well  justified. 

A  computer  program  was  implemented  using  the  VS  APL 
language  to  test  the  performance  of  the  LMS  adaptive  filter 
in  a  noisy  passive  environment.  The  array  size  used  in  the 
simulation  is  small  compared  to  most  modern  systems  but  the 
other  parameters  are  set  to  simulate  a  very  realistic  pas¬ 
sive  SONAR.  The  APL  language  is  used  because  of  its  large 
library  of  advanced  signal  processing  functions  and  its 
interactive  mode  of  operation  which  allows  for  rapid  program 
development.  The  major  disadvantage  of  the  APL  language  is 


its  slower  speed  of  execution  compared  to  running  a  compiled 
Fortran  program  that  performs  the  same  functions.  This  is 
not  a  serious  problem,  however,  for  the  research  of  this 
thesis . 

A .  PROBLEM  STATEMENT 

A  low  frequency  acoustic  signal  from  a  far-field  source 
is  received  by  a  planar  array  in  a  noisy  underwater  environ¬ 
ment.  The  noise  is  assumed  to  be  Gaussian  and  uncorrelated 
with  the  signal.  It  is  also  assumed  that  the  noise  between 
elements  is  uncorrelated.  The  planar  array  has  a  square 
structure  with  M  =  5  elements  in  the  x-direction  and  N  =  5 
elements  in  the  y-direction.  The  element  spacing  is  set  to 
one-half  the  wavelength  of  the  maximum  frequency  of  the 
system's  operating  range.  The  geometry  of  the  problem  is 
shown  in  Figure  17  [Ref.  18].  The  system  parameters  are: 

Signal:  A  cos(2-rrft),  where  A  is  the  amplitude  and 

f  is  the  frequency 

Integration  time:  0.5  seconds  (Tq) 

Frequency  resolution:  2  Hz  (1/TQ  =  fQ) 

Number  of  samples:  128  =  27  (L) 

Number  of  sensor  elements  in  the  x-direction:  5  (M) 

Number  of  sensor  elements  in  the  y-direction:  5  (N) 

Speed  of  sound:  1500  meters/second  (cQ) 

Sampling  rate:  256  samples/second  (f  ) 

s 

Frequency  range:  0-128  Hz 
Number  of  frequency  bins:  128 


Figure  17.  System  Geometry  {Ref.  18 :p.  19] 
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Elements  spacing:  7.5  meters  (d^  and  dy  equal) 

2 

Noise  model:  additive  white  Gaussian  ~  W{0,aN) 

Number  of  noise  samples  per  sample  function:  3200, 
i.e.,  L  *  m  x  N  =  128  x5  x5 

The  signal  used  in  the  simulation  is  a  100  Hz  sinusoid 
corrupted  by  white  Gaussian  noise.  A  total  of  128  samples 
are  taken  for  each  processing  period  of  0.5  seconds.  This 
corresponds  to  256  samples  per  second  which  satisfies  the 
Nyquist  sampling  theorem  [Refs.  12,16,17].  The  maximum 
observed  frequency  in  this  case  is  128  Hz.  The  100  Hz 
signal  will  center  in  bin  number  50  for  this  simulation. 
Frequency  bin  numbers  64-127  correspond  to  the  negative 
frequencies  [Ref.  16].  The  element  spacing  of  7.5  meters 
is  the  maximum  allowable  separation  for  the  100  Hz  signal 
to  avoid  grating  lobes  in  the  far-field  directivity  beam 
pattern.  The  speed  of  sound  is  the  speed  in  the  proximity 
of  the  planar  array. 

B.  SIMULATION 

Given  the  system  parameters  stated  in  Section  A,  and  a 
128  point  DFT,  a  100  Hz  sinusoidal  signal  will  be  centered  in 
frequency  bin  number  50.  The  logical  flow  graph  of  the  simu¬ 
lation  program  is  shown  in  Figure  18.  The  principle  of 
superposition  allows  different  frequency  bins  to  be  processed 
independently.  However,  if  the  same  frequency  is  emitted 
from  two  or  more  spatial  locations,  the  LMS  adaptive  filter 
will  lock  on  to  the  one  closest  to  boresight.  Complete 
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Figure  18.  Logical  Flow  Graph  of  the  Simulation 
Program  for  the  Low  Frequency  Passive 
Detection  Problem 
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documentation  of  the  simulation  program  for  this  case  can 
be  found  in  Appendix  B.  The  outputs  of  the  simulation  pro¬ 
gram  are  the  estimated  direction  cosines  u  and  v.  The  RMS 
error  is  defined  as: 


£ 


(AU2  + 


2  1/2 
Av  ) 


(3.1) 


where : 


Au  =  u  -  u 


(3.2) 


(3.3) 


where  u  and  v  are  the  actual  direction  cosines.  This  measure 
of  error  is  consistent  with  the  least-squares  criterion  used 
in  estimating  u  and  v.  If  the  estimate  of  the  spherical 

^  A 

coordinates  (9,0)  are  required,  the  following  transformations 

A  A  A  A 

will  transform  (u,v)  to  (9,0)  [Ref.  5]: 


9 (u, v) 


-1  ~  2 

sin  1  [  (up  + 


(vp] 


1/2 


(3.4) 


and 


AAA  i  A  A 

<p  (u , v)  =  tan'^Cv/u)  (3.5) 


It  should  be  noted  that  the  transformation  from  (u,v)  to 

A  A 

(9,0)  is  nonlinear.  A  particular  value  of  the  RMS  error 
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1 .  Orthogonal  Linear  Arrays 

This  configuration  uses  only  M+N-l  elements.  The 
parameters  used  for  the  simulation  are  as  follows : 

Signal  -  cos  2Tr(100)t  (bin  number  50) 

Number  of  iterations  (I)  -  75 

Initial  feedback  coefficient  (y  )  -  0.0005 

i 

Initial  feedback  coefficient  (y  )  -  0.0005 

y 

Scale  factor  for  y  and  y  (a  )  -  0.909  (see  Appendix  A) 

x  y  s  " 

Input  SNR  at  single  array  element,  i.e.,  SNR  at  FFT 
input  -  0  dB 

incident  angles  (9,cj>)  -  (55°, 35°)  (see  Figure  17) 

The  corresponding  direction  cosines  (u,v)  -  (0.67101,0.46985) 

Figures  19  and  20  show  the  convergence  characteristics  of 
the  complex  LMS  adaptive  filter  vs.  iteration  number  i.  The 
solid  horizontal  lines  are  the  true  values  of  direction 
cosines  u  and  v.  The  oscillations  in  the  beginning  of  the 
adaptation  are  due  to  the  large  initial  values  of  yx  and  y^,. 

The  feedback  coefficients  yx  and  y^  are  scaled  down  recur¬ 
sively  by  the  scale  factor  via  the  rule  described  in  Appendix 


A  which  is  rewritten  below: 
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Figure  19.  Convergence  Characteristics  of  the  LMS 

Adaptive  Filter  using  the  Orthogonal  Linear 
Arrays  Configuration  in  Estimating  the 
Direction  Cosine  u  at  an  Input  SNR  of  0  dB 
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Convergence  to  the  steady  state  occurred  in  fewer  than  35 
iterations  in  both  cases.  Figure  21  shows  the  error  in 
estimating  u  and  v  using  equation  (3.1)  vs.  iteration  number. 
The  use  of  a  constant  feedback  coefficient  requires  a  much 
longer  iteration  period  and  results  in  less  accuracy  (see 
Appendix  A) .  The  chosen  initial  values  of  ux  and  u  are 
outside  the  bound  of  the  convergence  coefficient  described 
in  Appendix  A.  However,  the  convergence  coefficients  are 
decreased  rapidly  using  equation  (3.10).  This  choice  of 
Ux  and  y  shows  that  the  LMS  adaptive  filter  with  decreasing 
convergence  coef f icient (s)  is  very  robust.  Table  6  summarizes 
the  simulation  results  for  a  particular  sample  function  of 
noise. 

2 .  Two-dimensional  Array 

This  version  uses  all  M  xN  elements  and  the  complex 
weights  are  not  assumed  to  be  separable.  The  parameters 
used  in  the  simulation  are: 

Signal  -  cos  2-n-(100)t  (bin  number  50) 

Number  of  iterations  (I)  -  75 

Initial  convergence  coefficient  -  0.001 

Scale  factor  (ag)  -  0.909 


Input  SNR  at  single  array  element,  i.e 
input  -  0  dB 


SNR  at  FFT 


REPRODUCED  AT  GOVERNMENT  EXPENSE 


ITERATION  NUMBER 


Figure  21.  The  RMS  Error  in  Estimating  u  and  v  Using 
the  Orthogonal  Linear  Arrays  Configuration 
at  an  Input  SNR  of  0  dB 


TABLE  6 


SUMMARY  OF  COMPLEX  LMS  ADAPTIVE  FILTER 
PERFORMANCE  FOR  THE  ORTHOGONAL  LINEAR 
ARRAY  CONFIGURATION  AT  SNR  =  0  dB 


Simulation  Values 


u 

0.671010 

V 

0.469846 

A 

u 

0.695675 

A 

V 

0.466038 

Au 

/\ 

=  u  -  u 

2.4665  xlo“2 

Av 

II 

<  > 

1 

< 

-3.808  xio'3 

e 

=  (Au2+Av2) 

0.024957 

0 

55° 

A 

6 

54.783° 

<D 

0 

in 

/V 

4> 

33.818° 

Incident  angle  (0,4))  -  (55°,  3  5°) 

Corresponding  (u,v)  -  (0.67101,0.46985) 

Figures  22  and  23  show  the  convergence  characteris 
tics  of  the  two-dimensional  array  configuration  of  the 
complex  LMS  adaptive  filter  vs.  iteration  number  i. 

Figure  24  plots  the  RMS  error  in  estimating  u  and 
v  (see  equation  (3.1)).  The  results  of  the  simulation  for 
a  particular  sample  function  of  noise  are  summarized  in 


Table  7. 
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TABLE  7 

SUMMARY  OF  COMPLEX  LMS  ADAPTIVE  FILTER  PERFORMANCE 
FOR  THE  TWO-DIMENSIONAL  ARRAY  AT  SNR  =  0  dB 


Simulation  Values 


0.67101 


0.46985 


0.68442 

0.46569 


1.341  x 10 


-4.16  xio 


1.404  xio 


55.876° 


34.232° 


3 .  Two-dimensional  Array  with  Separable  Weights 

This  configuration  assumes  that  the  complex  weights 
are  separable,  i.e.,  cd(m,n)  =  c(m)d(n).  The  simulation 
parameters  are  the  same  as  those  used  in  the  two-dimensional 
array  case  discussed  in  the  previous  section.  Figures  25 
and  26  show  the  convergence  characteristics  of  this  configura¬ 
tion.  Figure  27  shows  the  RMS  error  versus  number  of 
iterations.  The  summary  for  this  run  is  in  Table  8.  An 
alternate  way  to  illustrate  the  performance  of  the  complex 
LMS  adaptive  filter  is  a  plot  of  RMS  error  vs.  input  SNR  at 
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Figure  26.  Convergence  Characteristics  of  the  LMS 
Adaptive  Filter  for  the  Two-dimensional 
Array  with  Separable  Weights  in  Estimating 
the  Direction  Cosine  v  at  an  Input  SNR  of 
0  dB 
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TABLE  8 

SUMMARY  OF  COMPLEX  LMS  ADAPTIVE  FILTER  PERFORMANCE 
FOR  THE  TWO-DIMENSIONAL  ARRAY  WITH  SEPARABLE 
WEIGHTS  AT  SNR  =  0  dB 


U 

v 

A 

U 

*  /\ 
V 

Au 

Av 

£ 

e 

0 

A 

0 


Simulation  Values 
0.67101 
0.46985 
0.66304 
0.47353 
-7.97  x 10"3 
3.68  xio"3 
8.78  xio"3 
55° 

54.565° 

35° 

35.534° 


a  single  array  element  (SNR  at  FFT  input).  Figure  28  shows 
these  curves  for  all  three  configurations. 

D .  SUMMARY 

The  algorithm  that  assumes  separable  complex  weights 
shows  better  performance  over  the  prescribed  SNR  range.  The 
improvement  of  performance  of  all  three  algorithms  with 
increasing  SNR  is  evident  from  Figure  28.  Since  each  noise 
sample  function  has  a  total  of  3200  independent  samples,  the 
RMS  error  versus  SNR  curve  is  rather  smooth.  If  several  sample 
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REPRODUCED  AT  GOVERNMENT  EXPENSE 


Figure  28.  The  RMS  Error  in  Estimating  u  and  v 
Versus  Input  SNR. 

Curve  1:  Orthogonal  Linear  Arrays 
Configuration 

Curve  2:  Two-dimensional  Array 
Curve  3:  Two-dimensional  array  with 
Separable  Weights 
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functions  are  averaged,  the  slight  'humps'  in  Figure  28 
can  be  smoothed  out  further.  The  RMS  error  quantity  is 
analogous  to  that  of  a  'miss  distance'  on  a  rectangular  grid. 


IV.  PULSE  COMMUNICATION  [Ref.  18] 


The  simulation  results  of  two  cases  are  presented  in 
this  chapter.  They  are: 

Homogeneous  medium  case,  in  which  the  output  electri¬ 
cal  signal  data  at  each  array  element  is  produced  by 
the  ocean  communication  channel  simulation  computer 
program  [Ref.  19] .  Figure  29  shows  the  system  geometry 
of  this  case  [Ref.  18] .  Note  that  the  ray  path  from 
transmit  to  receive  array  is  a  straight  path. 

Inhomogeneous  medium  case,  in  which  the  ray  path  is 
bent  due  to  the  variable  sound-speed  profile.  Thus, 
the  apparent  direction  of  arrival  viewed  from  the 
receive  array  is  different  from  the  previous  case. 
Figure  30  shows  the  system  geometry  of  the  inhomo¬ 
geneous  medium  case  [Ref.  18] . 

From  the  analysis  of  Chapter  II,  it  can  be  seen  that  the 
LMS  adaptive  filter  should  be  able  to  phase  align  a  planar 
array  to  point  in  the  direction  of  arrival.  The  signals 
used  for  this  simulation  are  a  CW  pulse  and  a  LFM  pulse. 
These  are  very  common  signals  used  in  the  SONAR  environment. 
Information  is  carried  by  the  modulation  of  these  pulses. 

The  Fortran  program  used  to  simulate  this  problem  is  docu¬ 
mented  in  Appendix  C.  Blount  [Ref.  18]  studied  the  effect 
of  model-based  cophasing  on  the  probability  of  detection  of 
a  single  pulse.  The  amount  of  cophasing  is  determined  by 
the  system  geometry  and  deterministic  ray  bending.  It  was 
shown  that  by  applying  the  phase  weights  generated  by  con¬ 
sidering  those  factors,  the  performance  of  a  correlator 
receiver  was  improved  markedly.  Analysis  done  on  those 
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steering  phase  weights  showed  that  the  main  beam  of  the  array 
directivity  pattern  was  indeed  steered  to  the  direction  of 
actual  arrival  instead  of  line  of  sight. 

A.  TRANSMIT  WAVEFORMS  [Ref.  18] 

Two  types  of  waveforms  were  used  to  test  the  LMS  adap¬ 
tive  algorithm. 

1 .  Rectangular-envelope  CW  Pulse 

The  signal  presented  to  the  processor  is  a  quadra¬ 
ture  demodulated  complex  envelope  of  the  CW  pulse  [Refs.  5, 
18]  . 

K  jk(27rf  )t 

z(t)  *  J  z  e  °  (4.1) 

k— K  n 


where  denotes  complex  envelope. 

The  pulse  repetition  frequency  is  the  same  as  the 
fundamental  frequency  fQ  of  the  finite  (K  harmonics)  fre¬ 
quency  spectrum  from  which  the  pulse  is  ysnthesized.  The 
pulse  duty  cycle  is  arbitrarily  set  to  0.5.  The  complex 
Fourier  series  coefficients  zn  used  to  synthesize  the  com¬ 
plex  envelope  of  the  CW  pulse  are  obtained  from  a  closed-form 
expression  for  the  complex-valued  continuous  spectrum.  The 
Fourier  coefficients  are  obtained  by  evaluating  the  closed- 
form  expression  for  the  continuous  spectrum  at  discrete 
frequencies  corresponding  to  integer  multiples  of  the 


fundamental  frequency.  The  following  specific  transmit 
signal  parameters  were  used  in  all  CW  pulse  simulations 
Amplitude  (A):  40.0 

Duty  Cycle  (D) :  0.5 

Fundamental  Frequency  (f  ) :  200  Hz 

Harmonic  Values:  zq  =  20  exp  [j0°] 

z_^  =  zi  =  12.3324  exp  [j0°] 
z_2  =  z2  =  0*000  exP  I j 0 0 ] 
z_2  =  z3  =  4.244134  exp  [jl80°] 
z_4  =  z4  =  0.000  exp  [j0°] 
z_,_  =  z^  =  2.546479  exp  t  j  0  ° ] 

2.  Rectangular-envelope  LFM  Pulse  [Ref.  18] 

The  complex  Fourier  coefficients  used  to  synthesize 
the  LFM  pulse  are  found  using  a  procedure  similar  to  that 
used  for  the  CW  pulse  except  the  closed  form  expression  for 
the  complex-valued  continuous  spectrum  of  the  LFM  pulse  was 
found  by  using  the  method  of  stationary  phase.  Officer 
[Ref.  20]  describes  the  method  of  stationary  phase  as  does 
Papoulis  [Ref.  21]  who  also  provides  a  complete  description 
of  the  LFM  waveform.  The  following  transmit  signal  param¬ 
eters  were  used  in  all  LFM  pulse  simulations: 

Amplitude  (A):  40.0 

Duty  Cycle  (D) :  0.8 

Phase  Deviation  Constant  (B) :  2356.2  radians/volt 

Fundamental  Frequency  (f ) :  10  Hz 
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-  Harmonic  Values:  zQ  =  14.60593  exp  £ j 4 5 0 ] 

z_^  =  z\  ~  14.60593  exp  Ij21°] 

z_2  ~  z2  ~  ^4.60593  exp  [309°] 

z_3  =  z3  =  14.60593  exp  [jl89°] 

B .  PROBLEM  STATEMENT 

A  pulsed  signal  (CW  or  LFM)  is  sent  to  an  intended 
receiver  in  the  far-field.  The  signal  at  the  receive  array 
has  a  planar  wavefront.  The  direction  of  arrival  of  the 
incident  plane  wave  is  determined  by  applying  the  frequency 
domain  LMS  adaptive  filter  to  phase  align  (cophase)  signals 
at  all  sensor  elements  in  the  planar  array.  The  discrete 
time  signal  at  element  (m,n)  has  the  form: 

j^-(umdv+vndv) 

y  UTg,mdx,ndy)  =  zUTg)e  *  (4.2) 

where 


zUTg) 


K 


k=-K 


Zke 


j  k  2  it  f  o  t 


l :  time  index 

m:  element  index  in  the  x-direction 
n:  element  index  in  the  y-direction 
T:  sampling  period 


(4.3) 


u:  direction  cosine  with  respect  to  the  x-axis 


V 


direction  cosine  with  respect  to  the  y-axis 
d:  interelement  spacing  in  the  x-direction 

d:  interelement  spacing  in  the  y-direction 

z:  complex  Fourier  coefficient 
f:  fundamental  frequency 

K:  total  number  of  harmonics 

[Ref.  5]  shows  that  the  number  of  time  samples  needed  to 
completely  describe  equation  (4.3)  is  L,  where: 


L  >_ 

2K  +  1 

(4 

The  system  parameters  are 

• 

• 

CW 

LFM 

Integration  time  (Tq) : 

5  mS 

100  mS 

Frequency  resolution  (fQ) 

:  200  Hz 

10  Hz 

Number  of  samples  (L) : 

11 

7 

Number  of  sensors  (M) : 

5 

5 

Number  of  sensors  (N) : 

5 

5 

Sampling  rate  (f  ) : 

s 

2200  samples 

70  samples 

Number  of  frequency 

per  second 

per  second 

bins  (Q) : 

11 

7 

Element  spacing  d^,  d^: 

0.1229  meters 

0.1229  meters 

Carrier  frequency  (f  ) : 

w 

5  KHz 

5  KHz 

Noise  model:  Additive  white  Gaussian  noise  for 

both  cases 


Number  of  complex  noise 
samples/pulse : 


275 

97 


175 


The  carrier  frequency  f  is  assumed  to  be  known  a  priori.  It 
is  found  for  the  CW  pulse  that  the  best  direction  cosine 
estimates  are  obtained  by  processing  the  spectral  line  corres 
ponding  to  the  carrier  frequency.  This  is  reasonable  since 
the  signal-to-noise  ratio  at  that  bin  is  the  highest.  For 
the  LFM  pulse,  all  the  harmonic  lines  still  have  roughly  the 
same  magnitude  after  propagating  through  the  medium.  The 
accuracy  in  estimating  the  direction  cosines  is  about  the 
same  for  any  harmonic  line.  Therefore  the  spectral  line 
corresponding  to  the  carrier  is  processed. 

C .  RESULTS 

The  homogeneous  medium  case  is  considered  first  for  both 
CW  and  LFM  pulses,  followed  by  the  inhomogeneous  case. 

1 .  Homogeneous  Case 

The  parameters  for  the  system  geometry  (Figure  29) 

are: 

Speed  of  sound  (cQ)  :  1475  meters/second 

Depth  of  transmit  array  (YQ) :  1000  meters 

Depth  of  receive  array  (Yr) :  2500  meters 

Cross  range  (Xr  -XQ):  500  meters 

Line  of  sight  range  |r  - rQ | :  3000  meters 

True  spherical  angle  9:  31.81° 

True  spherical  angle  <J> :  -108.4° 

Direction  cosine  u:  -0.1666 


Direction  cosine  v:  -0.5000 


The  three  array  configurations  of  the  frequency  domain  LMS 
adaptive  filter  were  then  applied  to  the  CW  and  LFM  pulse 
cases  to  estimate  the  direction  cosines.  The  averaged  re¬ 
sults  for  one  pulse  corrupted  by  100  sample  functions  of 
noise  at  an  input  SNR  of  0  dB  at  a  single  array  element  for 
the  CW  pulse,  homogeneous  case  are  presented  in  Table  9. 

The  initial  convergence  coefficient  is  0.5  and  the  scale 

factor  a  is  0.909  for  both  CW  and  LFM  waveforms, 
s 

TABLE  9 

PERFORMANCE  OF  COMPLEX  LMS  ADAPTIVE  FILTER  FOR 
SPATIAL  RESOLUTION,  100  ITERATIONS,  INPUT 
-  SNR  =  0  dB  FOR  CW  PULSE,  HOMOGENEOUS  CASE 

Algorithm 

orthogonal  2-dimensional 

linear  arrays  with  separable  2-dimensional 


phase  weights 

array 

u 

-0.1666 

-0.1666 

-0.1666 

V 

-0.5000 

-0 . 5000 

-0.5000 

/V 

u 

-0.1714 

-0.1672 

-0.1670 

A 

V 

-0.4757 

-0.5003 

-0.5032 

Au 

-0.00048 

-0.0006 

-0.0004 

Av 

0.0243 

-0.0003 

-0.0032 

£ 

0.0247 

-0.0O0624 

0.0032 

9 

31.81° 

31.81° 

31.81° 

<p 

-108.4° 

-108.4° 

-108.4° 

9 

30.37° 

31.84° 

32.02° 

/V 

<P 

-109.8° 

-108.5° 

-108.4° 
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The  rms  error  versus  signal-to-noise  ratio  curves  for  all 
three  array  configurations  are  presented  in  Figure  31. 

The  simulation  results  for  the  LFM  pulse  in  the 
homogeneous  medium  at  an  input  SNR  of  0  dB  at  a  single  array 
element  are  summarized  in  Table  10.  A  total  of  100  differ¬ 
ent  noise  sample  functions  were  used  to  corrupt  the  received 
signal.  The  tabulated  values  in  Table  10  are  the  ensemble 
average  values.  Figure  32  illustrates  the  rms  error  versus 
SNR  plot  for  the  LFM  pulse  in  the  homogeneous  medium. 


TABLE  10 


PERFORMANCE  OF  COMPLEX  LMS  ADAPTIVE  FILTER  FOR 
SPATIAL  RESOLUTION,  100  ITERATIONS,  INPUT 
SNR  =  0  dB  FOR  LFM  PULSE,  HOMOGENEOUS  CASE 


Algorithm 


orthogonal 
linear  arrays 


2-d imensional 
with  separable 
phase  weights 


2-dimensional 

array 


u 

-0.1666 

-0.1666 

-0.1666 

V 

-0.5000 

-0.5000 

-0.5000 

A 

u 

-0.1096 

-0.1611 

-0.1462 

V 

-0.2695 

-0.4096 

-0.4079 

Au 

0.0570 

0.0055 

0.0204 

Av 

0.2305 

0.0904 

0.0921 

e 

0.2375 

0.0905 

0.0943 

0 

31.81° 

31.81° 

31.81° 

4» 

-108.4° 

-108.4° 

-108.4° 

A 

0 

16.91° 

26.11° 

25.68° 

-112.1° 

-111.5° 

-109.7° 

100 


RMS  ERROR 


Figure  31.  The  RMS  Error  in  Estimating  u  and  v  Versus 
Input  SNR  for  the  CW  Pulse/Homogeneous 
Medium  Case. 

Curve  1:  Orthogonal  Linear  Arrays 
Configuration 

Curve  2:  Two-dimensional  Array 
Curve  3:  Two-dimensional  Array  with 
Separable  Weights 
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SIGNAL  TO  NOISE  LEVEL  IN  DB 


Figure  32 


The  RMS  Error  in  Estimating  u  and  v  Versus 
Input  SNR  for  the  LFM  PUlse/Homogeneous 
Medium  Case. 

Curve  1:  Orthogonal  Linear  Arrays 
Configuration 

Curve  2:  Two-dimensional  Array 
Curve  3:  Two-dimensional  Array  with 
Separable  Weights 


The  performance  of  the  LFM  pulse  is  worse  than  that 
of  the  CW  pulse  since  the  number  of  samples  per  pulse  is 
seven  where  the  CW  pulse  has  eleven  time  samples  per  pulse. 
Consider  also  the  Fourier  coefficients  given  in  Section  A  of 
this  chapter.  The  magnitude  of  the  spectral  line  correspond¬ 
ing  to  the  carrier  for  the  CW  pulse  is  1.37  times  larger  than 
that  of  the  corresponding  spectral  line  in  the  LFM  pulse. 

However,  for  both  waveforms,  it  can  be  seen  from  Figures 
31  and  32  that  the  two-dimensional  array  with  separable 
weights  has  the  best  performance,  followed  by  the  two-dimensional 
array  and  orthogonal  linear  arrays,  respectively.  For  the 
CW  pulse,  accurate  spatial  localization  can  be  obtained  for 
SNR  greater  than  -3  dB.  For  the  LFM  pulse,  the  SNR  required 
is  about  3  dB. 

2.  Inhomogeneous  Case  [Ref.  18] 

The  parameters  for  the  system  geometry  (Figure  32) 

are : 

Speed  of  sound  {cQ) :  1475  meters  per  second 

Gradient:  0.017  per  second 

Depth  of  transmit  array  (y  ) :  1000  meters 

Depth  of  receive  array  (Y  ) :  2500  meters 

Cross  range  (X„  -X  ) :  500  meters 

r  o 

Line  of  sight  range  |r  - rQ  |  :  3000  meters 

True  spherical  angle  8:  30.10° 

True  spherical  angle  <p :  -109.4° 

Direction  cosine  u:  -0.1666 


Direction  cosine  v:  -0.4731 


Figure  30  shows  that  the  path  traveled  by  the  acoustic  rays 
is  bent.  The  LMS  adaptive  filter  is  able  to  resolve  the 
actual  direction  of  arrival  but  not  the  true  line  of  sight 
direction  to  the  transmit  array.  The  initial  convergence 
coefficient  u  is  set  equal  to  0.5  and  the  scale  factor  a 

s 

is  set  equal  to  0.909.  Both  the  CW  and  LFM  pulses  are  con¬ 
sidered  for  each  of  the  three  array  configurations.-  Table  11 
summarizes  the  performance  of  the  LMS  adaptive  filter  applied 
to  the  CW  pulse  case  at  an  input  SNR  of  0  dB  at  a  single 
array  element.  Figure  33  shows  the  decline  of  rms  error  as 
the  signal-to-noise  ratio  is  increased.  All  three  array 
configurations  are  included  in  the  plot  for  comparison. 

The  performance  for  the  LFM  pulse  case  is  tabulated 
in  Table  12  for  an  input  SNR  of  0  dB  at  a  single  array  ele¬ 
ment.  Figure  34  illustrates  the  performance  of  all  three 
array  configurations  versus  signal-to-noise  ratio  for  the 
LFM  pulse  case. 

D .  SUMMARY 

For  both  the  homogeneous  and  inhomogeneous  medium  cases, 
the  complex  LMS  adaptive  filter  performed  as  expected.  The 
array  configuration  that  assumed  separability  of  the  complex 
weights  consistently  demonstrated  better  performance  than 
that  of  the  orthogonal  linear  arrays  and  the  two-dimensional 
array.  The  superior  performance  can  be  attributed  to  the 
fact  that  equation  (2.24)  which  describes  the  reception  of 
a  plane  wave  by  a  planar  array  is  separable.  Therefore,  by 


TABLE  11 


PERFORMANCE  OF  COMPLEX  LMS  ADAPTIVE  FILTER  FOR 
SPATIAL  RESOLUTION,  100  ITERATIONS,  INPUT 
SNR  =  0  dB  FOR  CW  PULSE,  INHOMOGENEOUS  CASE 


orthogonal 
linear  arrays 

u  -0.1666 

v  -0.4731 

u  -0.1663 

v  -0.4527 

u  0.0003 

v  0.0204 

e  0.0204 

9  30.10° 

<0  -109.4° 

e  28.83° 


Algorithm 

2-dimensional 
with  separable 
phase  weights 

-0.1666 

-0.4731 

-0.1686 

-0.4782 

-0.0019' 

-0.0051 

0.0054 

30.10° 

-109.4° 

30.47° 


2-dimensional 

array 

-0.1666 

-0.4731 

-0.1651 

-0.4746 

0.0015 

-0.0015 

0.0022 

30.10° 

-109.4° 

30.17° 


-110.2° 


-109.4° 


-109.2° 
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Figure  33.  The  RMS  Error  in  Estimating  u  and  v  Versus 
Input  SNR  for  the  CW  Pulse/Inhomogeneous 
Medium  Case. 

Curve  1:  Orthogonal  Linear  Arrays 
Configuration 

Curve  2:  Two-dimensional  Array 
Curve  3:  Two-dimensional  Array  with 
Separable  Weights 
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TABLE  12 

PERFORMANCE  OF  COMPLEX  LMS  ADAPTIVE  FILTER  FOR 
SPATIAL  RESOLUTION,  100  ITERATIONS,  INPUT 
SNR  =  0  dB,  LFM  PULSE,  INHOMOGENEOUS  CASE 


Algorithm 


orthogonal 
linear  arrays 


2-dimensional 
with  separable 
phase  weights 


2-dimensional 

array 


-0.1666 


-0.1666 


-0.1666 


-0.4731 


-0.4731 


-0.4731 


-0.1035 


-0.1614 


-0.1417 


-0.2447 


-0.4058 


-0.3332 


0.0631 


0.0052 


0.0249 


0.2284 


0.0673 


0.1399 


0.2369 


0.0675 


0.1421 


30.10° 


30.10' 


30.10° 


-109.4° 


-109.4° 


-109.4° 


15.41° 


25.89° 


21.23° 


-112.93° 


-111.7' 


-113.0° 
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assuming  separable  complex  weights,  the  system  is  set  up  to 


match  the  physics  of  the  problem.  The  performance  of  this 


pulse  communication  system  can  be  enhanced  by  increasing  the 


pulse  width,  taking  more  time  samples  per  unit  time,  using 


high  resolution  spectrum  analysis,  and  enlarging  the  size 


of  the  array  by  adding  more  sensor  elements. 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  frequency  domain  LMS  adaptive  algorithm  has  been 
shown  to  perform  the  function  of  spatial  resolution.  The 
number  of  iterations  required  (approximately  35)  to  reach  a 
steady  state  is  found  to  be  much  fewer  than  that  of  a  com¬ 
parable  time  domain  adaptive  filter  [Refs.  2,4].  The  three 
modifications  made  to  the  original  complex  LMS  adaptive 
filter  [Refs.  1,4]  are: 

normalization  of  the  complex  weights  to  unity 
magnitude  after  each  update 

reduction  of  the  magnitude  of  the  convergence 
coefficient  for  each  iteration 

unwrapping  of  the  phase  weights 

It  has  been  shown  that  the  above  three  modifications  enable 

the  frequency  domain  LMS  adaptive  filter  to  be  applied  to 

a  multiple  element  array,  to  have  a  fast  convergence  rate 

and  robustness,  and  to  cover  the  entire  angular  range  of 

9  and  $  values. 

In  Chapter  III,  a  passive  low-frequency  signal  was 
generated  to  test  the  performance  of  the  frequency  domain 
LMS  adaptive  filter  in  the  presence  of  white  Gaussian  noise. 
The  low  frequency  problem  is  significant  in  the  area  of  long 
range  detection  and  possibly  long  range  control  and  communi¬ 
cation  for  ballistic  missile  submarines.  Some  form  of 
modulation,  whether  amplitude  or  angle,  will  have  to  be  used 


to  carry  the  information.  The  slow  data  rate  is  not  a  major 
drawback  since  only  certain  predefined  emergency  operation 
codes  will  be  transmitted.  The  simulation  results  indicate 
that  for  an  input  SNR  at  a  single  array  element  of  greater 
than  -3  dB,  accurate  bearing  and  depressive  angle  estimates 
can  be  obtained. 

In  Chapter  IV,  a  pulse  communication  problem  was  con¬ 
sidered.  Two  types  of  pulses,  rectangular-envelope  CW  and 
rectangular-envelope  LFM  pulses,  were  used.  The  possible 
applications  at  this  frequency  range  (5  KHz)  are  high  reso¬ 
lution  SONAR  imaging,  active  SONAR  detection,  and  communica- 
tion—between  submerged  vessels.  The  performance  of  the 
frequency  domain  LMS  adaptive  filter  was  again  tested  in  the 
presence  of  white  Gaussian  noise.  .  The  simulation  results 
indicate  that  accurate  bearing  and  depression  angle  esti¬ 
mates  can  be  obtained  for  input  SNRs  of  greater  than  0  dB. 
The  result  here  is  slightly  worse  than  that  of  the  passive 
low-frequency  case  due  to  the  shorter  data  length  used  for 
the  pulse  communication  waveforms.  For  the  simulations 
in  Chapter  IV,  the  number  of  time  samples  collected  is  the 
minimum  required  to  satisfy  the  Nyquist  sampling  theorem. 

The  following  three  configurations  of  the  frequency 
domain  LMS  adaptive  filter  were  considered: 
orthogonal  linear  arrays 

two-dimensional  array  with  separable  phase  weights 
two-dimensional  array 


In  all  simulations,  the  configuration  that  used  the  separa¬ 
bility  assumption  produced  the  best  results.  This  is  quite 
reasonable  since  this  configuration  uses  all  the  output  data 
from  the  available  elements  and  more  importantly,  its  assump¬ 
tion  of  separability  matches  the  physics  of  the  plane  wave 
signal  since  the  phases  of  a  plane  wave  in  the  orthogonal  x 
and  y  directions  are  separable. 

In  the  course  of  this  investigation,  some  interesting 

topics  for  further  research  revealed  themselves: 

exact  target  localization  (range,  depth,  bearing, 
and  depression  angle)  if  accurate  environmental 
data  can  be  obtained 

-  more  efficient  update  algorithm  for  the  convergence 
coefficient 

implementation  of  more  efficient  software 

application  of  the  steady  state  phase  weights  generated 
by  the  frequency  domain  LMS  adaptive  filter  to  Blount's 
[Ref.  18]  correlator  receiver 

addition  of  a  noise  reduction  system  before  spectral 
estimation 

modifications  to  make  the  system  jam-resistant 

applications  of  other  high  resolution  spectral  analysis 
techniques  to  produce  frequency  spectra  (such  as 
maximum  entropy  [Ref.  22],  autoregressive,  maximum 
likelihood,  etc.). 


APPENDIX  A 


THE  FREQUENCY  DOMAIN  LMS  ADAPTIVE  ALGORITHM 

A.  DERIVATION 

Consider  the  phase  aligning  system  shown  in  Figure  13. 
The  performance  objective  is  for  the  adaptive  filter  to  con¬ 
verge  to  a  set  of  phase  weights  such  that  the  array  output 
signal  will  match  a  reference  signal.  In  the  frequency 
domain,  the  output  of  a  linear  array  of  M  elements  at  iter¬ 
ation  i  can  be  written  as: 


Zi(q)  =  jf  l  ci(m)Y(q,m)  = 


m 


1  T 
M  ^(q) 


(A .  1 ) 


If  Z(q)  is  the  reference  signal,  then  the  error  between  the 
reference  and  the  array  output  is: 


‘i  =  Z(q)  -  Z^q) 


(A. 2; 


where  e^  is  a  complex  quantity. 
The  mean  square  error  is: 


eiei  = 


'ei 


=  [Z(q)  -Zi(q)  J  (Z(q)  -Zi(q)J*  (A. 3) 


Z  (q)  |  2  -Zi(q)Z*(q)  -Z*(q)Z(q)  +  1  Z±  (q)  |  2 


(A. 4) 
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I 

I 


(A. 9) 


— i+1 ^m)  =  +  2uiei^I  di (n)Y (q,m,n) ) * 

n 

— i+l (m)  =  -i(n)  +  2  iei{£  ci (m) Y(q,m,n) ) *  (A. 10) 

m 

where: 

ei  =  Z(q)  -  l  ci(m)  l  di (n) Y(q,m,n)  (A. 11) 

m  n 

Reference  4  studied  the  application  of  the  frequency  domain 

LMS  algorithm  to  a  split  array.  The  error  signal  is 

the  difference  between  the  outputs  of  the  two  arrays.  It 

is  possible  to  extend  this  idea  to  an  (M  xN)  planar  array 

by  generating  an  error  signal  e^Cn^n)  at  each  element  location 

where: 

ei(m,n)  =  Z(q)  -  cd^(m,n)Y(q,mfn)  (A. 12) 

for  all  m,n.  Although  this  scheme  is  not  studied  in  this 
thesis,  the  implementation  is  rather  straightforward.  One 
needs  only  to  replace  the  error  e^  with  e^Cn^n)  in  all  the 
equations  and  to  rewrite  the  code  to  generate  the  error. 

This  scheme,  however,  does  not  have  the  advantage  of 
averaging  noisy  data  points  from  all  of  the  elements  in  the 


array  to  come  up  with  an  estimate.  Therefore,  it  can  be 
assumed  that  the  noise  performance  will  probably  be  inferior. 

The  magnitude  of  the  reference  signal  Z (q)  is  generated 
by  averaging  all  (M  *N)  magnitude  spectra  whereas  the  phase 
of  Z(q)  is  taken  to  be  the  phase  of  the  reference  element. 

A  thresholding  is  done  to  determine  which  frequency  bins 
should  be  processed. 

B.  MODIFICATION  TO  THE  LMS  ADAPTIVE  ALGORITHM  TO  MAKE  IT 

USEFUL  FOR  SPATIAL  RESOLUTION 

Three  modifications  are  incorporated  in  the  basic  LMS 
adaptive  algorithm. 

1.  Normalization  of  Recursive  Complex  Weights 

The  components  of  the  recursive  complex  weights  (c^, 
d^  and  cd^)  must  be  normalized  to  have  magnitudes  equal  to 
one.  This  restriction  forces  the  algorithm  to  achieve 
minimum  error  by  adjusting  the  phase  weights  only.  The  esti¬ 
mation  of  direction  cosines  depends  on  the  linear  relationship 
of  the  phases  of  the  signal  across  the  array.  Figure  35 
shows  that,  without  the  normalization  performed  after  each 
complex  weight  update,  the  vector  sum  composing  the  estimate 
can  be  made  close  to  that  of  the  reference.  The  correct 
values  of  the  direction  cosines  can  only  be  achieved  if  each 
component  in  equation  (A.l)  has  the  same  phase  as  the  refer¬ 
ence  phasor.  All  the  vector  sums  shown  in  Figure  35  have 
the  same  resultant  phase  and  magnitude  as  the  reference  but 
the  phases  of  the  components  are  not  equal.  It  can  also  be 


Figure  35.  Different  Ways  of  Adding  Up  Vectors 
to  Obtain  the  Same  Resultant 
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seen  that  there  are  infinitely  many  ways  to  achieve  a  small 
error  without  the  normalization.  The  LMS  adaptive  filter 
implemented  for  this  thesis  has  a  unity  magnitude  across 
the  array.  It  is,  in  fact,  a  phase-only  filter.  With  the 
normalization,  each  component  is  forced  to  adapt  to  the 
phase  angle  of  the  reference.  This  is  the  only  way  that 
spatial  resolution  (based  on  the  amount  of  phasor  rotation 
of  adaptive  weights)  can  be  achieved. 

2 .  The  Convergence  Coefficient  (u) 

This  quantity  is  sometimes  referred  to  as  the  feed¬ 
back  coefficient,  signifying  the  analogy  to  control  systems 
[Ref.'  6].  For  nonstationary  environments,  a  constant  con¬ 
vergence  coefficient  is  usually  chosen  [Refs.  1,4,23]. 
However,  if  stationarity  can  be  assumed,  the  performance  of 
the  LMS  adaptive  filter  can  be  improved  by  using  a  mono- 
tonically  decreasing  scheme  for  implementing  the  convergence 
coefficient  [Ref.  11].  If  a  constant  convergence  coeffi¬ 
cient  is  used,  the  bound  on  the  choice  of  value  is  given 
by  [Refs.  12,24] : 


0  <  y  < 


max 


(A. 13) 


where 


I  I 

Amav  <  Tr [R]  =  l  R ( i , i )  =  £  R(0)  =  (I+1)R(0) 

max  i=0  i=0 


and  R  is  the  covariance  matrix  of  the  received  noise  cor¬ 
rupted  signal.  If  y  is  larger  than  the  maximum  determined 
by  equation  (A. 13),  oscillation  will  occur  [Ref.  25]  and 
the  LMS  adaptive  filter  will  not  converge.  If  y  is  too 
small,  then  the  rate  of  convergence  is  slow.  Figures  36,  37 
and  38  show  the  convergence  characteristics  for  different 
values  of  y . 

The  recursive  LMS  algorithm  is  based  on  the  stochas¬ 
tic  approximation  technique  [Ref.  11].  The  choice  of  y  is 
optimal  if  the  following  conditions  are  met  [Ref.  12] : 

y  i  >  0 

lim  y .  •  =  0 

i  -+oo  x 


r  2 

L  Ui  <  00 

i=l  1 

A  particular  choice  of  y^  that  meets  the  above  four  condi¬ 
tions  is  y ^  =  1/i  [Ref.  12],  The  effect  of  the  adaptation 
decreases  with  the  number  of  iterations  and  ceases  completely 
for  large  i.  Simulations  show  that  the  above  choice  of  y ^ 
is  better  than  using  a  constant  y,  however,  the  choice 


PERCENT  ERROR  IN  ESTIMATING  DIRECTION  COSINE 


CONVERGENCE  COEFFICIENT  TOO  LARGE 


Figure  36.  Convergence  Coefficient  (Constant  Value) 
Too  Large 


PERCENT  ERROR  IN  ESTIMATING  DIRECTION  COSINE 


PERCENT  ERROR  IN  ESTIMATING  DIRECTION  COSINE 


BEST  CONSTANT  CONVERGENCE  COEFFICIENT 


Figure  38.  Optimum  Convergence  Coefficient 
(Constant  Value) 
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where  0  <  ag  <  1  and  uQ  is  the  initial  value,  appears  to 
speed  the  rate  of  convergence  and  also  achieve  a  high 
degree  of  accuracy.  Equation  (A. 14)  does  not  meet  the 

CO 

condition  [  u.  =  »,  but  simulations  show  that  the  response 
i=l  1 

of  the  LMS  adaptive  filter  is  fast  and  accurate  for 

0.75  <  ag  <  0.9  5,.  Further  investigations  of  the  convergence 

characteristics  are  warranted.  The  current  implementation 

of  the  convergence  coefficient  results  from  experimentation 

with  various  values  of  a  to  achieve  the  best  estimates  of 

s 

direction  cosines.  With  this  scheme,  the  initial  value  of 
vu  can  be  set  quite  liberally  since  the  value  of  decreases 
geometrically . 

C.  PHASE  WRAP-AROUND  PROBLEM 

The  resolution  of  the  phase  wrap-around  problem  in  the 
adaptive  phase  weights  was  discussed  in  Chapter  II  exten¬ 
sively.  The  scheme  to  resolve  the  phase  ambiguity  depends 
on  the  proper  functioning  of  the  elements  adjacent  to  the 
reference  element.  In  the  event  that  some  of  those  adjacent 
elements  are  not  operational,  the  unwrapping  scheme  will 
have  degraded  performance.  However,  since  the  reference 
element  can  be  any  element  in  the  array,  it  is  possible  to 
shift  the  location  of  the  reference  element  to  a  region  of 
the  array  that  has  sufficient  adjacent  elements  that  are 
functioning  properly. 
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APPENDIX  B 

DESCRIPTION  OF  SIMULATION  PROGRAM 
FOR  THE  PASSIVE  DETECTION  CASE 

The  VS  APL  application  package  worksapce  (CMS  file 

ADAPTIVE  VSAPLWS }  contains  all  the  functions  necessary  to 

implement  the  simulation  discussed  in  Chapter  III.  The 

general  processing  flow  is  as  follows: 

generate  time  samples  of  a  plane  wave  signal  of 
frequency  f  incident  upon  a  (M  xN)  planar  array  with 
angles  of  incidence  (d,<p) 

add  white  Gaussian  noise  for  a  desired  SNR 

compute  the  discrete  Fourier  transform  (DFT)  of  each 
of  the  M  *n  time  sequences 

-  determine  the  spectral  line  with  the  largest  magnitude 
and  its  corresponding  frequency  bin  number 

apply  one  of  the  three  frequency  domain  LMS  adaptive 
filters  (orthogonal  linear  arrays,  two-dimensional 
array  with  separable  complex  weights,  or  two-dimensional 
array) 

compute  estimates  of  direction  cosines. 

Usage  of  the  functions  are  described  below. 

A.  SIG2D 

SIG2D  generates  planar  array  signal  at  each  element 
location  (equation  (2.28)): 

Syntax:  YN  -  AGl  SIG2D  AG 2 

YN  is  a  L  xM  xN  matrix  where  L  is  the  number  of  time 
samples,  M  is  the  number  of  elements  in  the  x-direction,  and 
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N  is  the  number  of  elements  in  the  y-direction.  AG1  is  the 
spherical  angle  9  and  AG2  is  the  spherical  angle  <p .  The 
speed  of  sound  (c) ,  the  signal  amplitude  (A) ,  the  frequency 
of  the  signal  (F) ,  the  interelement  spacing  (DX  and  DY) , 
the  sampling  interval  (T) ,  the  number  of  elements  (M,N) 
and  the  number  of  time  samples  (L)  can  be  changed  by  editing 
the  function. 

B .  NORRAND 

NORRAND  generates  independent  white  Gaussian  noise 
samples . 

Syntax:  NOISE  «-  K  NORRAND  N  Nl 

K  is  the  number  of  noise  samples  desired,  N  is  the  mean 
Nl  is  the  variance.  The  noise  array  NOISE  must  be  reshaped 
to  conform  to  the  shape  of  the  signal  generated  by  SIG2G 
by: 

NOISE  «-  L  M  N  p  NOISE 

The  standard  deviation  oN  =  1  is  necessary  to  scale  a  sample  func¬ 
tion  of  noise  with  zero  mean  and  variance  1  to  a  desired 

signal-to-noise  ratio.  The  signal  power  at  each  element  is 
2 

A  /2  and  the  noise  power  at  each  element  is  aN«  The  input 
signal-to-noise  ratio  at  a  single  array  element  is  then  given 
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cu 


(B.l) 


Solving  for  <?N  yields: 


A2  1/2 
°N  *  (2  (SNR)  ) 


(B.  2) 


where  A  is  the  amplitude  of  the  signal  and  SNR  is  the  numeri¬ 
cal  signal-to-noise  ratio.  Therefore,  a  noisy  signal  with 
the  desired  SNR  can  be  generated  as: 


RN  «-  YN  +  0N  X  NOISE  (B.3) 

where  NOISE  has  zero  mean  and  variance  1. 

C .  DFTWRT 

DFTWRT  computes  the  discrete  Fourier  transform  with 
respect  to  the  time  index  for  each  element  (equation  (2.42)). 


Syntax:  YK  «-  DFTWRT  RN 

YK  has  dimensions  L  xm  xN;  the  first  index  L  is  now  the 
total  number  of  frequency  bins.  RN  is  the  noisy  signal 
generated  by  adding  noise  to  the  output  of  SIG2D.  A  total 
of  (M  xN)  L-point  DFTs  are  computed  using  a  radix-2  FFT 
algorithm. 
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D .  ADPLMS 


ADPLMS  computes  estimates  of  the  direction  cosines  using 
the  frequency  domain  LMS  algorithm  for  the  orthogonal  array 
conf igurat ion . 


Syntax:  QTGT  ADPLMS  YK 

QTGT  is  the  frequency  bin  number  where  a  valid  signal  has 
been  identified  and  YK  is  the  output  of  the  function  DFTWRT. 

A  reference  signal  at  frequency  bin  QTGT  is  generated  by 
calling  the  function  REFGEN.  Direction  cosine  estimates 
u  and  v  are  computed  in  two  different  loops  since  in  general 
the  number  of  elements  M  is  not  necessarily  equal  to  N. 
Estimates  of  both  direction  cosines  are  generated  for  each 
iteration  under  the  names  UHAT  and  VHAT.  The  phase  unwrapping 
is  done  by  calling  the  function  DClDX  for  the  x-direction 
and  DC1DY  for  the  y-direction.  The  number  of  iterations 
(ITER)  and  the  initial  convergence  coefficients  (MUX,  MUY) 
can  be  changed  by  editing  the  function  ADPLMS.  The  scale 
factor  (SCMU)  for  decreasing  the  convergence  coefficients 
can  be  changed  in  the  workspace. 

E .  QFLMS 

QFLMS  computes  estimates  of  the  direction  cosines 
using  the  frequency  domain  LMS  algorithm  for  the  two-dimen¬ 
sional  planar  array  with  separable  complex  weights. 


Syntax:  QTGT  QFLMS  YK 


QTGT  is  the  frequency  bin  number  and  YK  is  the  output  of  the 
function  DFTWRT.  A  reference  signal  is  generated  by  the 
function  REFGEN.  Direction  cosine  estimates  u  and  v  are 
computed  every  iteration  of  the  LMS  adaptive  loop.  They 
are  stored  in  the  vectors  named  UHATQF  and  VHATQF.  The 
phase  unwrapping  is  done  by  using  DC1DX  and  DC1DY  in  the  x 
and  y-directions ,  respectively.  Recall  that  only  M+N  com¬ 
plex  weights  are  updated  for  this  configuration  because  of 
the  separability  assumption.  The  initial  convergence 
coefficient  (MUQ)  and  the  number  of  iterations  (ITER)  can 
be  changed  by  editing  the  function  QFLMS .  The  scale  factor 
is  named  SCMU  and  is  stored  in  the  variable  list  in  the 
workspace. 

F.  ADPLMS2D 

ADPLMS2D  computes  estimates  of  the  direction  cosines 
using  the  frequency  domain  LMS  adaptive  algorithm  for  a 
two-dimensional  planar  array. 

Syntax:  QTGT  ADPLMS2D  YK 

QTGT  and  YK  are  the  same  quantities  described  in  the  last 
two  functions.  A  reference  signal  at  frequency  bin  QTGT 
is  generated.  Direction  cosine  estimates  u  and  v  are  com¬ 
puted  for  each  iteration  and  stored  in  vectors  named. UHAT2D 
and  VHAT2D.  The  phase  unwrapping  is  accomplished  by  using 


the  function  DC2D.  The  number  of  iterations  (ITER)  and  the 
initial  convergence  coefficient  (MUG)  can  be  changed  by 
editing  the  function  ADPLMS2D .  The  scale  factor  can  be 
changed  by  assigning  a  different  value  to  SCMU  in  the  work¬ 
space.  The  use  of  a  different  unwrapping  function  is 
required  since  the  complex  weights  are  not  assumed  to  be 
separable. 

G .  REFGEN 

REFGEN  generates  a  reference  signal  at  a  particular 
frequency  bin  QTGT. 

Syntax:  CQREF  <-  QTGT  REFGEN  YK 

CQREF  is  the  reference  signal  used  in  the  frequency  domain 
LMS  adaptive  filter.  The  magnitude  of  CQREF  is  obtained  by 
averaging  the  magnitudes  of  all  (M  xN)  frequency  spectra 
in  the  frequency  bin  QTGT.  The  phase  of  CQREF  is  taken 
to  be  the  phase  of  the  reference  element  in  the  QTGT  fre¬ 
quency  bin. 

H.  DClDX 

DC1DX  unwraps  the  phase  weights  for  a  linear  array  in 
the  x-direction. 

Syntax:  UHAT  [ i ]  <-  N  DClDY  DV 
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UHAT [ i ]  is  the  direction  cosine  estimate  u  of  the  i 
iteration,  M  is  the  number  of  elements  in  the  x-direction 
and  CV  is  a  complex  vector  of  M  phase  weights  that  will  co¬ 
phase  the  incident  signal. 

I.  DC1DY 

DC1DY  unwraps  the  phase  weights  for  a  linear  array  in 
the  y-direction. 

Syntax:  VHAT  [ i ]  <-  N  DClDY  DV 

VHAT[i]  is  the  direction  cosine  estimate  v  of  the  i  iter¬ 
ation,  N  is  the  number  of  elements  in  the  y-direction  and 
CV  is  a  complex  vector  of  N  phase  weights  that  will  cophase 
the  incident  signal. 

J.  DC2D 

DC2D  unwraps  the  phase  weights  for  a  two-dimensional 
array. 


Syntax:  DC2D  CD 

CD  is  the  two-dimensional  phase  matrix  that  will  cophase 
the  incident  plane  wave  signal.  This  function  is  used  by 
the  function  ADPLMS2D. 

An  example  of  using  this  application  package  is  shown 


as  follows: 


YN  -  55  SIG2D  35 


NOISE  «-  3200  NORRAND  0  1 

NOISE  «-  128  5  5  p  NOISE 

RN  +-  YN  +  SCALE  NOISE 
YK  DFTWRT  RN 

QTGT  ADPLMS  YK 

QTGT  ADPLMS 2 D  YK 

QTGT  QFLMS  YK 


signal  generation 

noise  generation 

noise  generation 

noisy  signal 

transform  to  frequency 
domain 

estimate  direction 
cosines 

estimate  direction 
cosines 

estimate  direction 
cosines 


This  sequence  of  statements  generates  a  plane  wave  signal 
incident  upon  a  5  x  5  planar  array  with  angles  of  incidence 
0  =  55°  and  $  =  35°.  The  number  of  time  samples  for  each 
element  is  128.  A  noise  matrix  is  then  generated  and  added 
to  the  plane  wave  signal  and  the  discrete  Fourier  transform 
with  respect  to  time  is  taken  for  each  element  in  the  array. 
The  three  array  configurations  for  the  frequency  domain  LMS 
algorithm  are  then  used  sequentially  to  estimate  the 
direction  cosines  of  the  incident  plane  wave. 


APPENDIX  C 


DESCRIPTION  OF  SIMULATION  PROGRAM 
FOR  THE  PULSE  COMMUNICATION  CASE 

The  VS  FORTRAN  prgrams  were  written  to  implement  the 
pulse  communication  problem  discussed  in  Chapter  IV.  Two 
separate  programs  were  written  utilizing  essentially  the 
same  subprograms.  These  programs  handle  the  quadrature 
demodulated  complex  envelope  signals  generated  by  Vos' [Ref. 
19]  program.  The  programs  are  available  on  user  account 
0218P  at  the  Naval  Postgraduate  School,  Monterey,  California. 

A .  PROGRAM  ADBFM 

This  program  is  compiled  using  FORTVS  and  is  designed 
to  run  under  DISSPLA.  It  requires  a  storage  capacity  of 
1  Mbyte.  The  following  sequence  of  commands  should  be  used 
to  run  the  program. 

-  DEFINE  STORAGE  1  M 

-  I  CMS 

-  GLOBAL  TXTLIB  VALTLIB  VFORTLIB  CMSLIB  IMSLSP  NONIMSL 

-  GLOBAL  LOADLIB  VFLODLIB 

-  FILEDEF  04  DISK  fname  DATA 

(fname  is  the  filename  of  the  date  file) 

-  LOAD  ADBFM 

-  DISSPLA  ADBFM 

When  execution  begins,  the  user  will  be  prompted  to  enter 
the  desired  values  for  the  necessary  parameters.  These 


parameters  are  noise  status,  input  signal-to-noise  ratio  in 
dB  at  a  single  array  element,  number  of  iterations,  spec¬ 
tral  line  to  be  processed,  convergence  coefficient,  scale 
factor  for  the  convergence  coefficient,  and  the  choice  of 
one  of  the  three  array  configurations. 

For  each  array  configuration,  plots  of  the  estimates  of 
the  direction  cosines  are  generated.  The  plots  for  magni¬ 
tude  and  phase  of  the  difference  between  the  reference 
signal  and  the  estimate  are  also  generated. 

B.  PROGRAM  ERVSDB 

This  program  computes  the  rms  errors  for  various  input 
signal-to-noise  ratios.  If  a  plot  of  rms*  error  versus  SNR 
(dB)  is  not  required,  this  program  does  not  have  to  be  run 
under  DISSPLA.  The  following  sequence  should  be  used: 

-  DEFINE  STORAGE  1  M  (if  plot  is  required) 

-  I  CMS 

-  GLOBAL  TXTLIB  VALTLIB  VFORTLIB  CMSLIB  IMSLSP  NONIMSL 

-  FILEDEF  04  DISK  fname  DATA 

(fname  is  the  filename  of  the  data  file) 

-  LOAD  ERVSDB 

DISSPLA  ERVSDB  (if  plot  is  needed) 
or  START  *  (if  plot  is  not  needed) 

When  execution  begins,  the  program  will  prompt  the  user  to 

enter  an  initial  input  signal-to-noise  ratio  in  dB  at  a  single 

array  element.  It  will  then  ask  for  a  dB  step  size  such 

that  the  next  SNR  is  determined  by  the  current'  SNR  in  dB 

plus  the  dB  step  size.  A  total  of  nine  SNRs  are  allowed. 


For  example,  if  the  initial  SNR  is  -12  dB  and  the  step  size 
is  3  dB,  then  the  program  will  compute  rms  errors  for  the 
set  of  dB  levels  {-12,  -9,  -6,  -3,  0,  3,  6,  9,  12}.  The 
other  parameters  such  as  iteration  number,  initial  conver¬ 
gence  coefficient,  scale  factor  for  convergence  coefficient, 
and  the  spectral  line  to  be  processed  are  entered  when 
prompted.  The  program  will  then  ask  for  how  many  sample 
functions  of  signal  and  noise  are  to  be  averaged.  Simula¬ 
tion  results  show  that  the  average  of  50  to  100  sample  func¬ 
tions  are  sufficient  to  reduce  the  variance  of  the  direction 
cosine  estimates.  One  of  the  three  array  configurations  is 
then -chosen  by  the  user  to  estimate  the  direction  cosines. 

The  screen  output  of  this  program  is  ordered  pairs  of 
rms  errors  corresponding  to  a  particular  input  SNR  in  dB 
for  all  nine  specified  SNRs.  A  plot  of  rms  error  versus 
SNR  in  dB  can  be  generated  if  desired  (provided  that  the 
program  is  run  under  DISSPLA) . 
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